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To My Students 


Preface 


This book covers the finite element method and is intended to be used in courses 
taught in English or in bilingual courses. It is suitable for the knowledge system 
of undergraduate and postgraduate students majoring in engineering mechanics 
and similar majors. Most of the material in the book is based on the lecture 
notes of a bilingual course in computational mechanics given by the author. The 
lecture notes were used in undergraduate engineering mechanics majors from 
2017 to 2022 and the related graduate civil engineering, energy power, and 
mining engineering majors from 2019 to 2022 at China University of Mining and 
Technology (Beijing). The lecture notes have been revised and optimised several 
times to form this book. 

The book introduces the basic theory of the finite element method, formulates a 
finite element computation scheme for elasticity problems, and provides typical 
computer programs. To enable readers to quickly master the knowledge system and 
solution process of the finite element method, the book focuses on the basic con- 
cepts, algorithm theory, and key techniques of this method. From the perspective of 
reforming teaching methods and enriching the curriculum system, this book meets 
the following four requirements: 


(1) English instruction: The book is written in English. Simple and common 
English words are used to facilitate the interpretation of keywords, thus 
reducing the “double challenges” in the English barrier of non-native readers 
and the knowledge comprehension of new curriculum. 

(2) Systematic knowledge: The book only requires basic knowledge of elasticity. 
Starting with the simple weak form of the one-dimensional elasticity problem 
and the finite element solution process, it introduces the basic concepts and key 
techniques, namely, element, isoparametric transformation, numerical integra- 
tion, and computation scheme of two-dimensional and three-dimensional 
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VI Preface 


elasticity problems, to enable readers to quickly grasp the knowledge system of 
the finite element method. 

(3) Practice exercises: Carefully selected homework exercises at the end of each 
chapter and programming projects are well designed to form a set of exercise 
libraries matched with the knowledge topics of the book for testing the stu- 
dent’s understanding of the material and the corresponding exercise solutions 
are provided. Mastering the basic theory of the finite element method and 
having certain programming skills are the basic requirements of a computa- 
tional mechanics course. In addition to the theoretical knowledge, this book 
provides typical application programs so that readers can further practice by 
implementing the finite element method in specific problems. The electronic 
courseware related to this book can be available by contacting the author 
through the email: wangyl@cumtb.edu.cn. 

(4) Advanced applications: The book briefly introduces error estimation and the 
adaptive finite element method. This can serve as a reference for more advanced 
readers interested in high-performance computation and analysis techniques. 


This work was supported by the National Natural Science Foundation of China 
(Grant Nos. 41877275 and 51608301), Beijing Natural Science Foundation (Grant 
No. L212016), the Teaching Reform and Research Projects of Undergraduate 
Education of China University of Mining and Technology, Beijing (Grant Nos. 
J210613, J200709, and J190701), and the Yue Qi Young Scholar Project Foundation 
of China University of Mining and Technology, Beijing (Grant No. 2019QN14). 
Some of the material in this book refers to classical textbooks on the finite element 
method, and the pioneer authors and researchers are sincerely respected and 
appreciated. 

Because this book is restricted by the limited knowledge of its author, a few 
errors are unavoidable. The author hopes that experts, scholars, and other readers of 
this book will provide helpful suggestions for the book’s improvement. 


Dr. Yongliang Wang 

Department of Engineering Mechanics 

School of Mechanics and Civil Engineering 

State Key Laboratory of Coal Resources and Safe Mining 
China University of Mining and Technology (Beijing) 
July 2022 
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Chapter 1 


Introduction of Finite Element Method 


1.1 Development Process of Finite Element Method 


Understanding the behaviour and mechanisms of nature is fundamental in scientific 
research. Owing to the limited capacity of the human mind, we fail to grasp the 
complexity of the world in its entirety. Thus, engineers, scientists, or even economists 
subdivide a system under investigation into individual components, or ‘elements,’ 
the behaviour of which is easy to understand and reconstruct the original system 
from these components so that its behaviour can be naturally studied. 

In many cases, we can use a limited number of well-defined components to get 
an effective model, and such problems are called ‘discrete’. In other cases, the 
subdivision continues indefinitely, and the problem can only be defined by the 
fictitious concept of infinitesimals (which are not rigorously defined in standard 
mathematical analysis). This leads to differential equations or equivalent state- 
ments, which contain an infinite number of elements. Such systems are corre- 
spondingly called ‘continuous’. With the development of digital computers, 
discrete problems can usually be solved easily, even with a large number of ele- 
ments. Since the capacity of any computer is limited, continuous problems can only 
be solved accurately by mathematical manipulations. However, the mathematical 
techniques available for exact solutions are usually limited to oversimplified situ- 
ations. In order to solve realistic continuous problems, various discretisation 
methods have been put forward. All these methods involve an approximation that, 
ideally, as the number of discrete variables increases, it approaches the limit of the 
true continuum solution. 

Because most classical mathematical approximation methods and various direct 
approximation methods used in engineering belong to this category, it is difficult to 
determine the origin of the finite element method and the exact time of its invention. 
These approximate numerical methods definitely lay a crucial foundation for the 
generation of the finite element method, and promote the development of relevant 
numerical algorithms and numerical models in computational mechanics. The his- 
torical development of approximate numerical methods involved in the finite ele- 
ment method can also be obtained in some review papers. Figure 1.1 shows the 
evolution process, which leads to the concept of finite element analysis: 
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Trial functions 
Rayleigh, 1870 
Ritz,1908 


Weighted residuals Finite differences 
Gauss,1795 Richardson, 1910 
Calerkin, 1915 Liebman,1918 
Biezeno and Koch,1923| | Southwell,1946 


Variational methods 
Rayleigh,1870 
Ritz, 1908 


Structural analog substitution Piecewise continuous trial functions 
Hrenikoff,1941 Courant, 1943 
McHenry, 1943 Pragerand Synge, 1947 
Newmark,1949 Argyrisand Kelsey,1960 

Zienkiewicz and Cheung,1964 

Variational finite differences 

Direct continuum elements Varga,1962 
Turner et al.,1956 Wilkins,1964 
Clough.1960 Feng,1965 


Finite Element Method 


Fic. 1.1 — Historical development of approximate numerical methods involved in finite 


element method. 


(1) Mathematicians and engineers have different views on the discretization of 
continuous problems. Since 1870, some mathematicians have developed general 
techniques that can be directly applied to differential equations, such as the 
trial functions (Ritz, 1908; Rayleigh, 1870), variational methods (Ritz, 1908; 
Rayleigh, 1870), weighted residuals (Biezeno and Koch, 1923; Galerkin, 1915; 
Gauss, 1795). Furthermore, using the trial function in each partition domain, the 
piecewise continuous trial functions were developed (Zienkiewicz and Cheung, 
1964; Argyris and Kelsey, 1960; Prager and Synge, 1947; Courant, 1943). 

(2) The differential equations can be approximately expressed by difference equa- 
tions. In 1910, the classical difference method was established (Southwell, 1946; 
Liebman, 1918; Richardson, 1910). Then, based on the variational principle, the 
variational difference method was developed (Feng, 1965; Wilkins, 1964; Varga, 
1962). It should be emphasised that, in China, K. Feng of the Chinese Academy 
of Science also proposed a discretisation numerical method based on the vari- 
ational principle for solving elliptic partial differential equations in 1965. He is 
regarded as one of the pioneers in the basic theory of the finite element method: 
*Independently of parallel developments in the West, he (Feng) created a theory 
of the finite element method. He was instrumental in both the implementation 
of the method and the creation of its theoretical foundation using estimates in 
Sobolev spaces’ (Lax, 1993). 

(3) Some engineers and technicians directly used the method of analogy structure 
analysis to discretize the continuous domain, and the structural analog 
substitution was developed (Newmark, 1949; McHenry, 1943; Hrenikoff, 1941). 
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Subsequently, the method of direct continuum elements was developed, and the 
representative research works are from Turner et al. (1956) and Clough (1960). 
Turner et al. proved that a more direct but equally intuitive replacement of 
properties could be made significantly and more effectively by considering that 
small portions or ‘elements’ in a continuum behave in a simplified way. Clough 
used the term ‘finite element,’ which was born based on the view of ‘direct 
analogy’ in engineering, implying the direct use of a standard methodology 
applicable to discrete systems. This concept is extremely important, as it allows 
for improved understanding from the perspective of conceptual; meanwhile, it 
provides a unified method for various problems and the development of stan- 
dard computing programs. 


With the development of the above methods, various schemes for analysing 
problems of discrete nature have developed into a standard finite element method. 
When studying a structure, civil engineers first determine the force-displacement 
relations for each element of the structure and then assemble the structure by 
following a well-defined procedure to establish a local equilibrium at each ‘node’ or 
connecting point of the structure. The equations obtained can solve the unknown 
displacements. All of these analyses follow a standard pattern that is generally 
applicable to discrete systems. Thus, it is possible to define a standard discrete 
system, and this section focuses on establishing the procedures applicable to such 
systems. Using a unified treatment of standard discrete problems, the finite element 
process can be defined as a method of approximation to continuum problems as 
follows: 


(a) The continuum is divided into a finite number of elements whose behaviour is 
specified by a finite number of parameters. 

(b) As an assembly of its elements, the solution for the complete system precisely 
follows the same rules as those applied to standard discrete problems. 


The purpose of this book is to use the conventional and standard finite element 
method as a general discretisation procedure for constructing standard discrete 
problems from continuum problems. Furthermore, some high-performance finite 
element algorithms (such as the adaptive optimization algorithm) have been devel- 
oped; readers can further comprehend based on this basis of the standard method. 


1.2 Computation Procedure of Finite Element Method 


The reader should be aware that the finite element solution of the problem always 
follows the standard method, in which the process of solving any type of problem is 
always carried out by the following steps: 


(1) Define the problem to be solved by the governing differential equations. Based 
on the differential equations in the analysis of continuous systems, the equiva- 
lent integration form of the problem is constructed as virtual work, variational 
or weak form. 
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2) Select the type and order of the finite element. These elements and corre- 
sponding shape functions will be used in the analysis of discrete systems, which 
will be substituted into the equivalent integral form. 

3) Define a set of mesh for the problem. This involves the description of the 
distribution of nodes and elements and the description of the boundary 
conditions. 

4) Compute and assemble the matrixes from elements. The relationships in all 
elements are defined and assembled by the equivalent integral form of the 
problem; for example, the global stiffness matrix and load vector are set from 
each element using the element location vector. In this way, the continuous 
system is meshed and discretized to form an approximate discrete system. In 
this way, the differential equations representing continuous systems are trans- 
formed into algebraic equations of discrete systems. 

(5) Solve the system of algebraic equations. For widespread static problems, the 
algebraic equations are especially linear, which could be solved effectively by 
some well-developed mathematical techniques or specialized solvers. The solu- 
tions of the system of algebraic equations are solutions on all nodes. 

(6) Derive the solutions of variables in the global domain. Using the solutions on 

the nodes and the solutions of displacement, stress, and strain (derivatives of 

displacement) in each element domain can be obtained by further applying the 
interpolation of shape functions. 


Based on the above analysis steps, this book introduces the solution strategy and 
basic process of the finite element method. This book introduces the widely-used 
displacement-type finite element method. This method takes displacement u as the 
basic unknown variable. As shown in figure 1.2, in the continuous system, the gov- 
erning equations of mechanical behaviors are the differential equations Lu — 0 (where 
L is the differential operator) and the equivalent integral equation G (w, u) = 0 
(where G is the integral operator and w is an arbitrary function used) is obtained 
through equivalent transformation (such as weak form transformation). Then 
through discretization (using element and shape function), the algebraic equations 
Ku = f (where K and f are global stiffness matrix and load vector, respectively) on 
the approximate discrete system are obtained. Solving the algebraic equations can 


Continuous system Continuous system Discrete system 


Equivalent 
transformation 


Lu-0 D G (w,u)=0 D Ku-f 
L: Differential operator G: Integral operator K: Stiffness matrix 


Element and 
shape function 


Differential equations Discretization 


Integral equation Algebraic equations 


u: Unknown displacement w: Arbitrary function 


: Load vector 
Weak form J M 


u: Unknown displacement 


u: Unknown displacement 


Fic. 1.2 — Computation procedure from differential equations to integral equations and then 
to algebraic equations. 
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derive the approximate solutions of the basic unknown displacement u. This book 
provides, more specifically, the mathematical basis derived from these classic ideas. 


1.3 Main Contents of the Book 


This book introduces the basic theory of the finite element method for solving elastic 
mechanics problems and presents computer programs for finite element analysis. It 
is divided into nine chapters covering the above three aspects: elasticity mechanics 
theory, finite element methodology, and finite element program. In particular, the 
finite element methodology is the focus of this book, including basic techniques and 
computation schemes. Figure 1.3 presents a structural chart of the main contents 
and chapters: 


(1) In chapter 1, the historical development, as well as the computation procedure 
of the finite element method are briefly discussed, and a contents list is 
provided. 

(2) In chapter 2, the fundamentals of elasticity mechanics are introduced. 

(3) In chapter 3, the weak form of equivalent integration is presented. 

(4) In chapter 4, the elements and shape functions for one-, two-, and three- 
dimensional problems are discussed. 


Chapter 1. Introduction of finite element method 


Chapter 2. Fundamentals of Chapter 9. Programs of 
elasticity mechanics finite elementmethod 


Chapter 8. Error estimation 
and adaptive analysis 


Elasticity 


Chapter 3. Weak form of 
equivalent integration 


mechanics 
theory 


Chapter 7. Solutions of 
linear algebraic equations 


Chapter 4. Elements and ' 
ı Finite 

element 
| methodology 


shape functions 


Chapter 5. Isoparametric 
element and numerical 
integration 


Chapter 6. Finite element 
computation scheme of 


elasticity problems Finite 


element 
program 


a a a 


Basic techniques 


Fic. 1.3 — Structural chart of main contents and chapters in this book. 
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5) In chapter 5, the isoparametric element and numerical integration are 
introduced. 

6) In chapter 6, the finite element computation scheme of elasticity problems is 
described. 

7) In chapter 7, an effective numerical method for solving linear algebraic equa- 
tions is introduced. 

8) In chapter 8, the error estimation and adaptive analysis for high-performance 
algorithms are presented. 

9) In chapter 9, some typical computer programs for finite element analysis are 
provided, and their use is illustrated. 


1.4 Exercises 


1.4.1 Describe the treatment of ‘standard discrete problems’. 
1.4.2 Summarise the computation procedure of the finite element method. 


Chapter 2 
Fundamentals of Elasticity Mechanics 


This chapter covers equations for the elasticity theory and field theory, which are 
pretty common in diverse applications in engineering. First, the general 
three-dimensional elastic cases are studied and simplified to special two-dimensional 
problems. In the elastic equations given here, it is assumed that the strain is less 
than 1. It is assumed that the reader has got linear elasticity theory. Yet, for the sake 
of comprehensiveness, the basic equations for the different categories to be consid- 
ered are summarized. The reader can consult standard references and technical 
books about the theme for a more general discussion (Love, 2011; Ciarlet, 1988; 
Timoshenko and Goodier, 1970; Lekhnitskii, 1963; Sokolnikoff, 1956; Muskhelishvili, 
1953) and refer to the literature on large strain effects (Zienkiewicz et al., 2013a; 
Belytschko et al., 2000). In some situations, the indicial form of equations can be 
prior. Yet, in this part, the basic equations are expressed in the form of a matrix and 
later applied to establish finite element development. 

The basic equation of elasticity theory is explained by stress, strain, displace- 
ment, initial conditions, boundary conditions, and the constitutive relationship 
between stress and strain. First, each equation group of a general three-dimensional 
problem in a Cartesian coordinate system is designated, as displayed in figure 2.1. 

In addition, two-dimensional forms are considered. The issues of the two 
dimensions are as follows: 


(a) Plane stress case. There is only nonzero stress in the problem plane here and 
no stress in the direction orthogonal to the thin plate, as exhibited in 
figure 2.2a. 

(b) Plane strain case. It is assumed that the strain perpendicular to the plane 
under consideration is zero. This may occur in a prism, as shown in figure 2.2b, 
where the load perpendicular to the plane remains unchanged. 

(c) Aaxisymmetric case. In the cylindrical coordinate system r — z — 0, the angle 
0 in the plane considered is constant, as displayed in figure 2.2c. It is 
assumed that all components of stress, strain, and displacement are only related 
to r and z. 
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x(u) 


Fic. 2.1 — Three-dimensional elasticity problems. 


x r(u) 


(a) Plane stress (b) Plane strain (c) Axisymmetry 


Fic. 2.2 — Two-dimensional elasticity problems. 


2.1 Displacements 
In terms of the three-dimensional issue, the field of displacement is provided by: 


u (2, y, 2) 
u(x) = 4 v(zx,y,z) p, (2.1) 
w (£, y, z) 


where the position is given in Cartesian coordinates xz, y, z, which may also be 
expressed in the form of a matrix as: 


X= yp. (2.2) 
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For the case of two-dimensional plane strain and stress, the displacement field is 


illustrated by: 
_ Jus y) 
ud) eph (2.3) 


and only the first two components of equation (2.2) are required. In the end, in the 
axisymmetric case, the r — z plane is defined by cylindrical coordinates, thus 


r 
x= { : \ (2.4) 
with the displacements shown by: 
_ J ulrz) 
u(x) = { vifa) \ (2.5) 


The difference between the latter two is the coordinates that are applied: z and 
y represent Cartesian coordinates, while r and z represent cylindrical coordinates. 


2.2 Strains 


There are six independent strain components in a three-dimensional problem. They 
are arranged in order and expressed in the form of a matrix, that is: 


E= | €z Ey êz Vay Yyz Ya]. (2.6) 


where [.]" is the transpose of the matrix. The form is called Voigt notation in the 
mechanical literature (Zienkiewicz et aL, 2013b). This is an approach to writing 
symmetric second-order tensors with a group of simplified components. Strain is 
symmetrical, where Yay = Vyss 7, = Yap and Ya = Vrz; hence, the nine components 
are declined to six by using the Voigt notation mentioned above. Here, in terms of 
the problems of two dimensions considered, the last two components are forever 
zero. Hence, four components of e are supposed to be taken into consideration, which 
will be mentioned in the following section. 


2.3 Stresses 


In a body with three dimensions, there are nine stress sectors acting on a volume 
element, as exhibited in figure 2.3. The components Tap Tyr, Tyz, Tey Taz, and Tyg are 
known as shear stresses, while o,, c,, and o, are known as normal stresses, the shear 
components are symmetrical (balanced by angular momentum or moment on the 
cube), hence 


Tay = Tye, Ty = Tey and o Tu = tu. (2.7) 
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Fic. 2.3 — Stresses in three-dimensional solid analysis. 


With the same situation as strain, six components are applied to convey stress, 
which is ordered and expressed in the form of a matrix, given by: 
T 
06 —|os Oy Oz Try Ty Tx] 


(2.8) 


In the problems of the plane with two dimensions, only four stress components 
are taken into consideration, shown as: 


o= |r Oy o; Bel. (2.9) 


In the problems of plane stress, a priori is known as o; — 0. 
In terms of the issue of axial symmetry, the stress components are defined, as 
displayed in figure 2.4, as follows: 


6—[o, 90, 09 £s]. (2.10) 


Fic. 2.4 — Stresses in axisymmetric solid analysis. 
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2.4 Geometric Equations 


The strains of small deformation problems are calculated by displacement and can 
be conveyed in the form of a matrix as: 


e = Su, (2.11) 


where § denotes the differential operator matrix, and u represents the field of 
displacement. In a problem of three dimensions, the strain-displacement relation is 
expressed as: 


ð 

az 0 0 

0 m 0 
br oy 
by 0 d " 

s=] c SHE (2.12) 

Jay ð 0 0 " 
yz Oy Ox 
Vex 

ages 

Oz Oy 
ð à 
az ° Oe 


In order to give thought to the three kinds of two-dimensional issues in a coin- 
cident way, four strain components are included in e, hence, 


0 
Ox 
Er O 0 
0 I 
E He Oy HL 0 = S,u +8; (2.13) 
Ez U £z 
. 0 0 ^ 
i a 8 
Oy Ox 


in terms of the problems of the plane (where e; denotes that plane strain is zero, note 
that it is not suitable for plane stress), and 


— 


8 
5; 0 
" 5.18 
M mM OF MET (2.14) 
€0 v 
= 0 
n T 
dope 
Oz Or 
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in terms of the issues of axisymmetry, as displayed in figure 2.4. The difference 
between the types of two-dimensional problems only exists in plane stress problems 
g, and the existence of component e; in the case of axisymmetry. 


2.5 Constitutive Equations 


All the above equations have nothing to do with the material properties under 
consideration. In order to implement this theory, the relationship between stress and 
strain should be established, which may reflect the properties of different materials. 
The relations are called stress-strain equations or constitutive equations, which can 
arise in complex forms, yet in this procedure, the simplest linear elastic case is taken 
into consideration. 

Under this assumption, the equations of stress-strain for a linear elastic material 
are given as: 


o = De (2.15a) 
or 
e= D'o, (2.15b) 


where the matrix D is called the elasticity matrix, and its inverse, D'!, is the 
elasticity matrix of compliances (Timoshenko and Goodier, 1970). 

In terms of materials of elasticity, the stress can also be inferred from the energy 
storage function W displayed for the strains, therefore 


_aw 


E (2.16) 


It is applicable to both linear and nonlinear materials, hence, it can be utilized to 
generalise the above forms. In terms of the linear elastic equations in equa- 
tion (2.15a), the stored energy function is as: 


W(e) = 5eDe. (2.17) 


The complementary energy U(e) can also be defined from which the strain can be 
derived as: 


QU 
=— 2.18 
e= (2.18) 
where in the linear elastic case, 
1 
U(e) = 39 Dio. (2.19) 


The application of a complementary or stored energy form also guarantees that 
D and D™ are symmetric. 
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In isotropic materials, there is no change in the matrix of elasticity when a 
coordinate alternation is used. The general expression of isotropic materials is 
written with six stress and strain terms. In terms of isotropic materials, any two 
independent elastic constants can be used. 

For instance, under Cartesian coordinates, the form is offered by: 


Er 1 -v —v 0 0 0 Or 
Ey —v 1 —v 0 0 0 Cy 
e | 1|-v -v 1 0 0 0 6; 
yf E|0 0 0 2ü0-v) 0 0 tw d (220) 
Vaz 0 0 0 0 2(14 v) 0 Tiz 
Ys 0 0 0 0 0 2(14 v) Tur 
'The elasticity matrix of moduli is obtained by inverting: 
(1— v) v U 0 0 0 
v (1— v) v 0 0 0 
OE v v (1— v) 0 0 0 
nec] i 0 0  (1-20/2 0 0 y (221) 
0 0 0 0 (1— 2v)/2 0 
0 0 0 0 0 (1 — 2v)/2 


where d = (1 + v)(1 — 2v). In order to keep all the parameters of materials positive, 
the form of d limits the assigned material parameter values of v. Hence, 


1 
-l<v< z (2.22) 
The value range of the parameter is - 1 < v < 1/2. 
In terms of isotropic materials, the form of a two-dimensional problem is con- 
stituted using four stress and strain terms via ignoring the last two columns and 
rows in equation (2.20). Via Cartesian coordinates, the form turns. 


Ex 1 -v —v 0 Or 
& | 1|-v 1 -v 0 Oy 
& (| El|l—-v —v 1 0 o, ( (2.23) 
ds 0 0 0 21+) ] (ta 


In the example of plane stress, c, must be set to zero to calculate the matrix 
D. Thus 


& = = 2 6). (2.24) 


The inverse of equation (2.23) can be obtained. 


Or 1 UV 0 0 Er 
Oy Z E v 1 0 0 Ey 
o (7 ü-45j|0o00 0 B; 2) 
Tay 0 0 0 (1-v)2 Yay 
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The third row and column of array D are filled to obtain the stresses. 

In the plane strain and axisymmetric problems, the inverse is directly obtained 
via equation (2.23), as all stress components are nonzero. The following expression 
can be derived in Cartesian coordinates: 


Or (1— v) U U 0 £y 
oy | FE v (1— v) U 0 £y 
06; ( d v v (1— v) 0 £; [" (ean) 
Try 0 0 0 (1 = 2v)/2 Yay 


which is consistent with the problem of three dimensions if the last two rows and 
columns of each array in equation (2.22) are ignored. It is noticed that the case 
where e, is zero is handled by only substituting this value when the stresses are 
computed, yet c; is always nonzero unless the Poisson’s ratio v is zero. 


2.6 Equilibrium Equations 


2.6.1 Three-Dimensional Problems 


'The linear momentum or equilibrium equations for the solid behaviour of three 
dimensions can be described in Cartesian coordinates, given by 


Oo, | Ot; | Ota 


dx | Oy Oz a" 
ny Ooy Oty 
a + +b, = 2.27 
Ox Oy Oz je 22) 
Ot; Oty | Oo, u 
Ox Oy * Oz imi 


The stress is symmetric for angular momentum (moment equilibrium). Using 
only the independent components of the symmetric stress offered in equation (2.8), 
the equilibrium equations, which are in Cartesian coordinates, can be expressed by: 


STo+b = 0, (2.28) 


where, in the system of Cartesian coordinates, S denotes the differential operator, 
which is identical to that for strain in equation (2.11), and b represents the vector of 
body forces, shown by: 


b—[b, b, b]. (2.29) 


2.6.2 Two-Dimensional Plane Stress and Strain Problems 


In the two-dimensional plane problems, t,,, Tz», and b, are ignored. It is noted that o; 
does not relate to the z coordinate. Hence, the stress matrix is expressed by equa- 
tion (2.9), and the two equilibrium equations left can be provided as: 
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Oo, $ Oty; ipei 
Ox Oy 
(2.30) 
Ot, , Wy ap ecd 
Or ^ Oy T 
By noting that t,, = t,,, the equilibrium equations can be displayed by: 
T — 
S o +b = 0, (2.31) 


where S, is identical to that in equation (2.13), and b and u have only two 
components. 


2.6.3 Two-Dimensional Azisymmetric Problems 


In axisymmetric problems, the stress is displaced in equation (2.10), and body force 
is displayed by b = |b, Gs Via the stresses on the volume element expressed in 
figure 2.3 and summing the forces in the directions of r and z, the two equilibrium 
equations are shown as: 


Oo, _or— 90 | OT or 
T 


+ b, = 0 
Or T Oz 
P (2.32) 
Ór Oz AUT 


By moment equilibrium (angular momentum), z,, = T; is obtained again. Note 
that in the axisymmetric problem, the equilibrium differential operator is shown as: 


Lo 0 1 o 
sre Or r r OZ 


a ð 9 
Oz Or 


(2.33) 


and is unequal to 8T. This distinction is owing to the system of curvilinear 
coordinates. This is a disadvantage if the differential equation form is used, yet for 
approaches applied in the finite element analysis, it is exhibited that the distinction 
vanishes, and expressions in the curvilinear coordinate system allow significant 
simplifications. 


2.7 Boundary Conditions 


The equilibrium and strain-displacement equations are valid for all body domain 
points, denoted by Q. Boundary conditions (denoted as T) are appointed for each 
point on the solid surface. T'here are two kinds of boundary conditions to be thought 
of: (a) displacement boundary conditions and (b) traction boundary conditions. The 
boundaries can be sorted into two classifications, containing the displacement 
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conditions on the domain I, and the traction conditions on the domain I';. Hence, the 
entire boundary is the collection of the displacement domain and traction domain. 

'The condition of the displacement boundary is appointed at each point of the 
boundary T, as: 


u—u(x), xcI,, (2.34) 


where @ is the given value and x is the boundary point. 
The traction boundary conditions are designated for all boundary points on T; 
and provided for the following stress equation: 


t=G'e=i(x), xeT, (2.35) 
where, as to three-dimensional problems, G” denotes the matrix provided by: 


Cale du 0 ma n% 0J, (2.36) 


0 0 n O n, n 


with n, n, and n; being the direction cosines for the outward-pointing normal n to 
the boundary I. Note that matrices G and S are with identical nonzero structures. 
In two-dimensional plane problems, G is reduced to: 


T Ny 0 0 Ny 
Gi = E ni 2 (2.37) 


where n, and n, are the components of the unit outward normal vectors on the 
boundary. In the case of axial symmetry, n, — n, and m, -— m, are obtained, 
meanwhile, G, is also nonzero-like G,,. 

Both types of boundary conditions are vectors, constraining all of these points on 
the boundaries. It is promising to designate some components by displacement 
conditions and other components by traction conditions. 


2.8 Exercises 


2.8.1 Summarise the basic variables and equations of three- and two-dimensional 
elasticity. 

2.8.2 Using the formulation presented in this chapter, provide the basic variables 
and equations of axial deformation for the one-dimensional elastic cantilever beam 
shown in figure 2.5. The coordinates of the two ends are z — a and z — b, and the 
axial uniform load is b+ 


Fic. 2.5 — Axial deformation of one-dimensional elastic cantilever beam. 


Chapter 3 
Weak Form of Equivalent Integration 


3.1 Weak Form of Equivalent Integration for Differential 
Equations 


The previous chapter lists governing partial differential equations of elastic mechanics 
and the boundary conditions. The general solutions obtained by this method are 
rarely suitable and even challenging to solve in science and engineering. Therefore, we 
need to find other solutions; for example, we can use other forms to characterize the 
identical problem in order to get a more accurate approximate solution with this 
method. In this chapter, we will introduce the weak form of differential equations. 
Through the following steps, we can obtain a weak form for a set of differential 
equations: 


(1) Multiply each differential equation by an appropriate arbitrary function. 
(2) Integrate the product over the space domain of the problem. 

(3) Reduce the order of the involved derivatives using the integration by parts. 
(4) Introduce the boundary conditions reasonably. 


In this chapter, we only integrate into the space domain, therefore it is unnec- 
essary to express the time dependency and initial conditions. Its integral time 
domain can also be considered (Katona and Zienkiewicz, 1985; Wood, 1984; 
Zienkiewicz et al., 1984), but it is not necessary to study here. It can construct a 
weak form when decomposing any differential equation by using the above steps. An 
arbitrary function is one that can take any conceivable value. The arbitrary function 
can use a polynomial, a trigonometric function, a Heaviside function, a Dirac delta 
function, or any other function to take any conceivable value. The characteristic of 
arbitrary function is virtual displacement, which may be the reader already knows. 
Here, let’s consider a weak form for the one-dimensional equilibrium equation of 
elasticity first, which can help us understand the above process. 


3.2 Weak Form of One-Dimensional Elasticity Problems 


First, a one-dimensional elasticity problem with a single component z in a Cartesian 
coordinate system is considered a representative example. The displacement 
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u(x) can be expressed as a function of the coordinate z. It can be assumed that the 
problem is defined in a domain Q, and also have some interval that a < x « b. So we 
write the displacement to u(z) and found that all strain components to be zero as 
defined in the previous chapter except 


Ou 
Eg E 3.1 
Jz (3.1) 
Reduce the system of multiple equilibrium equations into a single equation: 
00, 
b, — 0, 3.2 
a to (3.2) 
and derive the single constitutive equation as follows: 
Or = Ey, (3.3) 


where E is Young's modulus of elasticity. 
'The following boundary conditions will be introduced: 


u-"uü or i,—í,-m,0, onzc- a,b, (3.4) 


where n, is the unit outward normal, which is +1 at b and —-1 at a. 

We can combine the above equations by substituting the constitutive and strain— 
displacement equations into the equilibrium equations. Then we can get the equi- 
librium equation expressed by displacement as: 


ð ðu 
E (522) NT (3.5) 


where E can be a function of z. There is only a single dependent variable u(x), it is an 
irreducible form of the differential equation. 

The set of equations (3.1)-(3.4) is referred to as the strong form of the 
one-dimensional elasticity problem. We can use these equations to formulate a set of 
discrete equations obtained by the finite difference method (Hildebrand, 1987; 
Mitchell and Griffiths, 1980; Forsythe and Wasow, 1960). However, a better way is 
to use the steps described above to introduce an integral weak form of the equations. 
Then the discrete equations are derived based on this weak form. 

First, in domain Q described by an interval a < x « b, we introduce an arbitrary 
function w(x). By multiplying the equilibrium equation (3.2) by this function, we get 
that 


g(@, u, Fx) = w(z) (v =) =0. (3.6) 


It is more convenient to have a form equal to zero because if w(x) is zero 
everywhere except where z = ze, then for equation (3.6) to be zero, at z = z,, it has 
to satisfy the differential equation. By repeating all points, it is a conclusion that at 
all points in the domain differential, the equation is satisfied, provided that w(x) is 
arbitrariness. Then, in the next step, we obtain an integral form that is also zero by 
integrating over the domain: 
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CGi E f "T IC + Flan =0. (3.7) 


In the third step, we can integrate the stress term into parts as follows: 


Oo, un z=b _ ðw 
] 5 dz = w(z)o;] |i. nz 


z 
= [uw(z)o;]|;-s — [w(z)os]l;-a — n 


ow (3.8) 
= [w(z)e,n;]|,-s + [w(2) ox Nz] t-a - [ Zoar 
Q Ox 
= [w(z)o, ny] r- f oda 
Q Ox 
where I' is the boundary of Q, and we have n,(a) = — 1 and n,(b) = 1. The 
boundary term can be expressed as follows by tractive force: 
[w(2)o2Me]|p= [w(z)t]]r- (3.9) 


We also denote T, by the part of the boundary where u = u and T; by the other 
part of the boundary where t, = tz; and the entire boundaries are I = l',UI;. With 
this symbol, we can express the weak form of the equilibrium: 


G(w, u, Or) = n watas- [ oe da wt;|r, + wte|p 
Q [e] Ox : 


ðw z 
= | w(2)b,dz J aade + wti;|r, + wt; 
ncc E: lr, + wills, 
0 


(3.10) 


For simplicity, by the special properties of the arbitrary function, specify the 
place of displacement, we set w|ļr, = 0 to eliminate the need of t, 


G(w, u, or) = f w(2)beae— | oda wte|p = 0. (3.11) 
Q o Ox i 


To obtain a pure displacement or ‘irreducible’ form of the equations, we can 
introduce the constitutive equation and strain displacement equation. Thus, we got 
that 


G(w, u) = [wo oraz- f es EP aa wtz|r,= 0. (3.12) 


'This form is known as a displacement model or an irreducible form. By com- 
paring equation (3.12) with equation (3.5), we found that the strong form requires 
the second derivatives of u and w with respect to the coordinate xz, but the weak form 
only needs the specification of the first derivatives. Construction of approximate 
solution further needs to be based on weak form. 
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3.3 Finite Element Computation Based on Weak Form 


3.3.1 Galerkin Method 


We multiply a specified set of functions by unknown parameters to express the 
displacement u(z, t), thus the approximate solution is constructed. Similarly, we use 
the same amount of specified functions and arbitrary parameters to write an arbi- 
trary function. These can be expressed as follows: 


N 
u(x) = U(x) = V ^ o, (2)as + uz) 

pa (3.13) 
w(x) = w(x) = 5 Wml) Wm 


We assume that in the above form, both qo,(x) and w,,(x) are zero at the 
specified displacement boundaries. Then it specifies the function u;(z) as any 
function satisfying the displacement boundary condition. For instance, if the dis- 
placement boundary condition is u(L) = d on the domain 0 < x < L, the function 
can be defined as: 


The Galerkin method is considerably older than the finite element method. The 
latter mainly uses locally based (element) functions in the expansion of equa- 
tion (3.12), but the general procedures are identical. Since this process produces 
equations with integral form, these equations can be obtained by weighted sum- 
mation of the subdomains; all weighted residual approximations are collectively 
called generalised finite element method. Occasionally, we find that it may be useful 
to use ‘local’ and ‘global’ trial functions. 

The approximate Galerkin solution to the elasticity problem can be obtained by 
inserting equation (3.13) into equation (3.12): 


D= Yun [6.6 dz — Yu, f 2. 


m=1 m=i 


N do, do; — 
A à; ig 


dx | an 
(3.14) 


+ 5 WinW au(x Jt, zr, = m 0 


Because the functions, Pp, Ym, and pz are known, the integrals can be evaluated 
as: 


dY m do, 
o dr dz qs, 


d 
fon [ varie~ [is Vm p 5 ddz + we lr, 
, en ae 


Kn = 


(3.15) 
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The array Ky, is the stiffness matrix, and fm is the specified load matrix. Because 
the parameters Wm are arbitrary, the expression multiplied by each wm must be zero. 
'This leads to the following set of equations: 


N 
S Kmnan = fm, m — 1,2,..., N. (3.16) 
n=1 


The problem is thus expressed in the matrix form of the stiffness equation: 


Ka =f, (3.17) 
where 
Ky Ky > Kın fi 
T a Kaz s "€ f 
Km Km © Kyn fy 


The formal solution to the problem is now given by: 
a—K f. (3.18) 


However, the inverse of K is never computed in practical situations. Instead, the 
solution can be obtained either by directly solving equation (3.17) or using an 
iterative approach. Once the parameters a, are known, we can use equation (3.13) 
to compute the approximation for the displacement, and that for stresses from: 


N 
_ _ do, des(z) 5 
o, (r)—-E am zb» dz an + da d]. (3.19) 


Example 3.1. Solutions of Galerkin method for a one-dimensional elastic cantilever 
beam. 


For instance, we take a static problem with a length of 10 units and E = 1000, as 
shown in figure 3.1. 


t=25 


Fic. 3.1 — One-dimensional elastic cantilever beam in example 3.1. 
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The loading is given by: 


p 10 for0<x<5 
"* |0 for5<x<10’ 


and the boundary conditions are 
u(0)— 0 and £1,(10)- 25. 


The weak form of this problem is 


10 5 
| PU 900 ode — n w(z)10dx — w(10)25 = 0. (3.20) 
0 Ox Ox 0 


We consider an approximate solution given by: 
N ENa N gNm 
ü(x) = » (55) a, and w(x) = /» (|) Wm. 
Note that u;(z) = 0 for this problem. The solution to the weak form is given by: 
N 10 m-—1 n-1 5 m 
E x x 
10 (=) (=) dz| an = I (=) 10dz+ 25, 3.21 
E mnis id x} a MEST z4 (3.21) 


where m — 1, 2,..., N. The stiffness matrix and the load vector are 


10 m- n—2 T 
deum n l0mn(—) dz = QUUM. 


$ gm 100 Jy 7! 
m= —] 1 +25 = + 25. 
f f Go) eee sia) x 


For example, if M — N — 2, the stiffness matrix and the load matrix are 


— [Ak ka! 100 100 | [Ah f 75/2 
ME D. e 7 E na £A hal ` (9:22) 
Using equation (3.18), the unknown parameters can be obtained as follows: 
|. fa _ f 0.62500 
iE { a j = { —0.25000 \ ne) 


The unknown parameters obtained from the above solutions are substituted into 
equation (3.13), and the displacement function can be obtained as follows: 
2 x n x xz 2 
xila) =S (2) a = Z x 0.625 (=) —0.25 
u(x) & ü(x) P» (5) On = 36 * + 0 x ( ) 


= —0.00252? + 0.06252 


(3.24) 
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Using equation (3.18), the stress function is obtained as follows: 


par « dos = —5z + 62.5. (3.25) 


Or 4 dr ”" 


Similarly, the solutions using one to five terms are presented in table 3.1, in which 
the symbol N represents the number of terms. 


TAB. 3.1 — Parameters for one-dimensional elastic cantilever beam in example 3.1. 


N a a a3 a4 di 

1 0.37500 

2 0.62500 —0.25000 

3 0.78125 —0.71875 0.31250 

4 0.78125 —0.71875 0.31250 0.00000 

5 0.73437 —0.25000 —1.09275 1.64063 0.65625 


We can find that there are some critical aspects in constructing approximate 
solutions from the results of the above problem. First, we note that the solutions 
using three and four terms are the same. This means that in the weak form, terms 
that contribute to the solution can be discarded. Another conclusion is when the 
solution is constructed in a global manner, it does not result in convergent behaviour 
for individual parameters. Thus, it is necessary to use more stable approximation 
procedures, such as the finite element method (this will be discussed later). The 
solutions using one to three terms for displacement and stress are obtained and 
shown in figure 3.2. It can be observed that the displacement at the free end is the 
same, no matter how many terms are used. This is only more common and valid in 


x 
(a) (b) 


Fic. 3.2 — Displacement and stress solutions in example 3.1 based on Galerkin method using 
N terms: (a) displacement and (b) stress. Notes: The symbol N represents the number of 
terms, and the symbol Exact represents the exact analytical solutions. 
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one-dimensional static problems and less applicable in higher-dimensional problems. 
The displacement solution converges faster than the stress solution, but we found 
that some points are more accurate than others. These points are called supercon- 
vergent points, which play a significant role in future research on error estimation 
and adaptive mesh refinement. 


3.3.2 Finite Element Computation 


It is a more convenient way to construct the approximating functions y, and q, by 
dividing the domain to be analyzed into small regular-shaped regions. For example, 
between a and b, we can divide a one-dimensional region into ‘M’ finite small 
segments by defining a set of N points x; such that 


Tj— 0, mj«zrj,j1 and zy - b. 


For a one-dimensional problem, it is also possible to have each increment define a 
finite element domain (or, more simply, an element), and groups of points define 
nodes. The element division and node distribution are the basic components of the 
finite element method for describing what is referred to as the finite element mesh or 
simply the mesh for the problem. 

By using the subdivision above, we can define a simple group of continuous 
polynomial approximating functions: 


0, X-—Ta 
T — X1 
— — —, _ G15 X S Gi 
Ti — G1 
9; — MEM (3.26) 
Tii T 
—— —— —, QTL Tig 
Tip] — Ti 
0, T Tii 


'The plot of these functions and their first derivative is shown in figure 3.3. Only 
when the function is continuous in z, but the first derivative is piecewise continuous, 
and the discontinuities are located on the nodes we call it Co function. In subsequent 
chapters, we considered selecting appropriate C$ functions for higher spatial 
dimensions. Here, we just study the description of general forms in one dimension. 

It can be seen from the figure the end functions can be used as the spatial form of 
the @; functions if necessary. 

We assume that y; = @,, thus, over the domain, all integrals defined in the weak 
form functional can be obtained. Indeed, by noting that, we can usually evaluate the 
integrals on each element individually 


[vee 2 3 [AM Ode 308 (odas (3.27) 


Considering any interval |z; z;,1] shown in figure 3.4, we note that it is defined 
by the same two local functions Mı and N5, which are called the shape functions for 
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$, h EMEN NS ha (1-99 
$, E ii. EEEE E dp /dx I-A 
Qı MES A NEEE T R 52» 2 oe 
$., NN NOE do. /dx e——À à 
Py RET d db Jie — —— ss) 
A - X WEE Xei Xs Hp dp o Ap 0 0 eem Xi Xn 
(a) (b) 
Fic. 3.3 — One-dimensional finite element approximation for q,;: (a) functions and 


(b) derivatives. 


(a) 


Fic. 3.4 — One-dimensional finite element shape functions: (a) functions and (b) derivatives. 


the element. We also define the local node coordinates in each element as xf and 2; 


to simplify the notation. 


With the above notation, the displacements and the arbitrary weight function 


also can be expressed as: 


(3.28) 


In each element, we define the local coordinate system as y = x — xf, and let 
he = zj — xf be the element length; the shape functions can be expressed as: 


x 


—1-—— 
he 


and they are the same for all elements 

given by the following formula: 
dN dM . 
dz dz |— 


Ni (a7) 


he 


and M(x) = (3.29) 


he’ 


. The derivatives of the shape functions are 


dN, dN, 1 


ane ae dal he” 


(3.30) 
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Thus, we can write the approximation of the weak form as: 


G( tv, à) = G,(à, à) — Gil, à) — ü(x)tulr,, (3.31) 


where the contribution to internal stress is Gi and the contribution of body loading 
to force is Gr. Using equation (3.28), each term is defined as follows: 


d Ni 
M h ^ 
> A A ae A dz’ dN, nd { 
G, (i, à) = Y ^ [ae ag E, da/4 ^L 4, 
(tv, à) - [ài as) f dN» | dz dz’ S Us 
da’ 


(à, à) = Yl às] n i { "i Joa (3.32) 


Each element can be evaluated as: 


dN, 
K= [ da! g [s emus - 1 E 
0 dM dz da’ Ky, Ky, (3 33) 
da’ 


- fO as = (411 
r= [Gp 


For the shape functions given in equation (3.29), if we make E, and b, in each 
element constant, the matrices are 


Afi A "HEINE 


3.4 Global Assembly from One-Dimensional Elements 


In the one-dimensional beam model shown in figure 3.5, N — 1 elements (‘M small 
finite segments) and N nodes are used. 


Element number @ Q 


Node number 1 2 3 N-1 N 


Fic. 3.5 — Element and node numbers of one-dimensional problem. 
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Then, we define the relation between the local and global node numbers for each 
element, as indicated in table 3.2. 


TAB. 3.2 — Local-to-global node numbering for two-end elements for one-dimensional elastic 
beam. 


Element number 1 2 3 Ut N-1 
Local node number 1 1 2 3 Ut N-1 
2 2 3 4 Ue N 


According to the node numbers of each element, the corresponding element 
location vector can be provided to determine the relation between local and global 
node locations as follows: 


Ael 2}7, AZ eIqg9 3)" e[N-2i N}T. (3.35) 
To globally assemble one-dimensional elements, the node numbers are shown on 


the left and upper sides of the matrix, and the expanded form of the stiffness matrix 
is given through the element location vector: 


Node 1 2 3 TH N-—1 N 
number 1 ri Ki Y A ù 
Ka Ko LAB 2 0 Ug 
3 0 Ka (Ki + Ki) Kj (is 
t 1 Ky — (Kx +i) Kaj | twa 
Ky Kə UN 
fi 
b 
"E 
- : (3.36) 
fua 
Iv 


Similarly, the node numbers are shown on the left of the load vector, and the 
expanded load is given by: 


Node 1 f a P 
number 2 2 Ti 
3 3 Tt 
` = eee N (3.37) 
N-1 ma N-2 fN 
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The summation indicated in equation (3.36) and subsequent expressions result 
in a standard linear problem, with the final stiffness equations being 


Ku=f (3.38) 
where 
Ky Ky +: Kin fi 
Kı Kə Kin hb 
K=] . i and f= 
Kn, Kyo © Kyn fn 


Now, the formal solution to the problem is: 


u—K f. (3.39) 


But it is difficult for us to compute the inverse of K in practical situations. 
Instead, the solution is obtained either by directly solving equation (3.38) or by an 
iterative method. Once the displacements at the nodes u are known, the approxi- 
mation of the displacement for each element may be computed as: 


(a!) = Ny (a) uo + Ny (a) E, (3.40) 


and that for stresses from as: 


ô (x)= E 


D(L) —,(dN(z),., dNX(z),. 
a -s( as UR. (3.41) 


3.5 Treatments on Boundary Conditions 


This section is concerned with the boundary conditions of the one-dimensional beam 
model. Especially, we discuss the displacement and traction boundary conditions for 
representative one-dimensional elastic problems. Other issues can be treated simi- 
larly at the corresponding nodes. 


(1) Displacement boundary conditions 


Using the above finite element form, because the parameters are now all physical 
values, it may be easy to implement the displacement boundary conditions. That is, 
they satisfy 


alta) = Ta. (3.42) 


Thus, to impose the displacement condition u(0) = Ñ, we set à; = u and rewrite 
equation (3.36) as: 
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1 0 0 uy ui 
0 Kx + Kj Ki . 0 [Ur h- Kat 
0 Kj, (Ka + Kj) Ki w ( "E 
Ky ^ (Kt + Ky) KG | | êv- Ivi 
KA KA UN fn 


(3.43) 


This is equivalent to setting tw; = 0. To simplify, the row and column of known 
displacement at the corresponding nodes can be deleted, thus reducing the size of 
the stiffness matrix. 


(2) Traction boundary conditions 
Similarly at x = L, impose the traction condition only requires the modification: 


These physical characteristics highlight the obvious advantages of using the 
finite element method to approximate variables. 


Example 3.2. Solutions of finite element method for one-dimensional elastic can- 
tilever beam. 


The example shown in figure 3.1 is considered again here. The domain was 
divided into four equal elements, as shown in figure 3.6. 


Fic. 3.6 — Four-element mesh for one-dimensional elastic cantilever beam in example 3.2. 


For example, using four elements, the stiffness matrix, and the load matrix are 


Be} 1 —1 1000| 1 —1 400 —400 
1 2 3 4 _ E a = 
diario ei dd s E 1 | 2.5 E 1 | E 400 | 


1 1 1 1 12.5 
1 2 5 z = 
f =f = shah b= zx iox25{ b= Sizs 
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pape iad ib=Zxox2s{ tey (3.45) 


Using the element location vector, the global stiffness matrix may be derived as 
follows: 


Node 1 2 3 4 5 
number 

1 [Kh Ki 

2 | Ka (KatKì) Kh 

3 Ka (KS, + Kj) : Kj 

4 Kj, (Ke + Kh) Kj 

5 Ki Ko 

400 —400 
—400 400+ 400 —400 (3.46) 


= —400 400 + 400 —400 
—400 400+ 400 —400 
—400 400 
400  —400 
—400 800  —400 
= —400 800  —400 
—400 800  —400 
—400 400 


Similarly, the load matrix may be derived as follows: 


Eb 1 (f fh 12.5 12.5 
2 |$ 14 pe 12.54+12.5 25 
3 d frase Ptfps=2 12540 $=? 125 (3.47) 
4 A z +f 0+0 0 
5 (5 fi 0 0 


The matrix form of the stiffness equation can be expressed as follows: 


400 — —400 ity 12.5 
—400 800  —400 ity 25 
—400 800  —400 ig 9 —4 12.5 (3.48) 
—400 800  —400| | à 0 
—400 400 | | as 0 


Then, the displacement boundary conditions à; = wu, = 0 can be introduced into 
the stiffness equation, and the global stiffness equation can be rewritten as follows: 


1 0 ity 0 
0 800 —400 în 25 
—400 800  —400 ig $= % 12.5 (3.49) 
—400 800 —400| | à 0 


—400 400 Us 0 
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To simplify, the row and column of known displacement at the corresponding 
node 1 can be deleted, thus reducing the dimension of the stiffness matrix: 


800  —400 tg 25 
—400 800  —400 ûs (| J 12.5 
—400 800  —400 ia ( — 0 (a50) 
—400 400 Us 0 


The traction boundary conditions f, — f; + ¢,(10) = 25 can be introduced into 
the load equation, and the global stiffness equation is then rewritten as follows: 


800  —400 tg 25 
—400 800  —400 ûs | _ J 12.5 
—400 800 —400 ta f — 0 Gan) 
—400 400 Us 25 


Using equation (3.39), the unknown displacements at the nodes are obtained as 
follows: 


fin 0.15625 
ts (J 0.25 
iy (~~ ) 0.3125 (uus) 
tis 0.375 


The unknown displacements on the nodes obtained from the above solutions are 
then substituted into equation (3.40), and the displacement function for each ele- 
ment can be obtained as follows: 


/ 


/ 
Elementl:  à(z)- (: e =<) x0+ = x 0.15625 = 0.06252" (3.53a) 


i / 
Element2: @2(z’) = (1 = x) x 0.15625 + 7 x 0.25 = 0.03752’ +0.15625 (3.53b) 


€ 


: zd 
Element3:  à?(z)-— (1 — =<) x 0.25 + 7” 0.3125 = 0.0252’ + 0.25 (3.53c) 


/ 


/ 
Element4:  à*(3)— (: E x) x 0.3125 4- = x 0.375 = 0.0252/ + 0.3125 (3.53d) 


e € 


Using equation (3.41), the stress function for each element can be obtained as 
follows: 
Ou! (x^) 
Ox! 


1 1 
Element1: ól(z)-—E = z( yx 0+ V 0.15625 = 62.5 (3.54a) 


€ 
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mad R i 1 1 
Hemed Sa ta a( = X 0.15625 + — x 025) = 37.5 (3.54b) 


€ € 


1 1 
" ^35. = =s = 
Element3: G°(2') = am z( XX 0.25 4- PRA 0.3125) 25 (3.54c) 


€ € 


Q^ (x 1 1 
fia: Pia e( = X 0.3125 + > x 0.375) =25  (3.54d) 


The solution above is shown in figure 3.7, in which the symbol n represents the 
number of finite elements. Similarly, it is also shown the solution uses eight elements. 
Comparing the two solutions, we note that the displacements converge more quickly 
than the stresses once again, but they have some superconvergent points. 


x 
(a) (b) 


Fic. 3.7 — Displacement and stress solutions in example 3.2 based on finite element method 
using n elements: (a) displacement and (b) stress. Notes: The symbol n represents the number 
of finite elements, and the symbol Exact represents the exact analytical solutions. 


3.6 Exercises 


3.6.1 Describe the implementation steps for the weak form of equivalent integration. 
3.6.2 Using the Galerkin method, for example 3.1 with E(x) = Ex, b,(x) = byg, 
compute the global stiffness matrix K and load matrix f. 

3.6.3 Using the finite element method, for example 3.2 with E(x) = Ex, b,(x) = b,z, 
compute the global stiffness matrix K and load matrix f. 

3.6.4 Summarise the functions of the element location vector. 


Chapter 4 
Elements and Shape Functions 


To compute the element matrices, it is necessary to select appropriate shape func- 
tions. Previously, we introduced the basic strategy of finite element computation 
based on the weak form using a one-dimensional element with two nodes, the sim- 
plest low-order linear element. In this chapter, we introduce a high-order Lagrange 
element for one-dimensional problems. Further, for two-dimensional problems, 
the triangle and rectangle elements are provided; for three-dimensional problems, 
the tetrahedron and hexahedron elements are also presented. 


4.1 One-Dimensional Lagrange Element 


4.1.1 Linear Element with Two Nodes 


This linear element has only two end nodes corresponding to the endpoints of a line 
segment, as shown in figure 4.1, and the displacements at the nodes are &àf and (5. 
Such elements are called one-dimensional elements with two nodes. 


Element e 


Lr— ÀÀÁ n4 


Nodenumber 1 2 


Fic. 4.1 — One-dimensional element with two nodes. 


For this one-dimensional problem, we write the displacement as: 


S 

Q 

i: 

I 

Me»: 

E 
ul 

II 


N (a!) uo + Ny (a^) àE. (4.1) 
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Previously, a local coordinate system in each element was defined as 2' = x — af, 
and the element length is he = x} — xj; the shape functions are given by: 


N(a)-1-L and N(x’) =— (4.2) 
he he 
and are the same for every element. The derivatives of the shape functions are 
dN, dN, 1 dN, dN, 1 
== — d o = = —. 4. 
de dr he dx dv OW (4.3) 


In equation (4.1), the displacement function is linear in z; therefore, the element 
is called a one-dimensional linear element. 


4.1.2 Higher-Order Lagrange Element 


It is difficult to describe complex displacement changes using the aforementioned 
linear elements. Accordingly, higher-order shape functions and elements with more 
nodes are introduced. A higher-order element with more nodes is shown in 
figure 4.2. 


Element e 


ie ü d 


ü ue 
e o o o 2 


Node number 1 3 4 


Fic. 4.2 — One-dimensional higher-order element. 


The displacement is written as: 


u* a = Y NACE) üt = Ni (E) WU + M(E) à+ + NCE) às, -1<E<1, (4.4) 


a=1 


where n is the total number of given functions used in the element, à; are unknown 
parameters to be determined, and € is a new local coordinate. The local coordinate 
selected here is for the convenience of numerical integration in the local coordinate 
system in the next computation procedure. The local coordinate C is related to 
the global coordinate x as follows: 


= 3 Nal) of = M(E HNE H N(z-lxéxi (45) 
a=1 


The global and the local coordinate in the above transformation are in 
one-to-one correspondence, which is not necessarily linear. Only when the shape 
functions are linear, as in the case of the two-node elements obtained previously, is 
the transformation linear. 
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To ensure that the unknown parameters to be determined are the displacements 
at the nodes, the following equation should hold true: 


1 = 
N49) = 8a 0 ues (4.6) 


where č, is the local coordinate with position a. Then, through the values of 
equation (4.4) on the node coordinates, we can get that the unknown parameters are 
exactly the displacements on the nodes. 

In order that equation (4.6) can be satisfied, the following simple construction 
for higher-order shape functions may be used, based on the Lagrange interpolation 
formula: 


T (6 =e) 

i(S) = ee 
AL, (ča = Čo) 
bza 


£ (é — EE — &)C-J(£— £y (E — Eagar ME En) (4.7) 


(oy = E1) (Ea m &)C i (Ca E Sant) (ha m Gat: d nes m £s) ' 
where the order of the polynomial P(E) is p = n-1. The coordinates of internal 
points of €, are uniformly located between the endpoints. 
For one-dimensional elements, the shape functions are defined as follows: 


NalS) = (8) (4.8) 


Further, the completeness condition requires that w°(€) should contain any 
constant c (displacement of rigid body), yielding. 


Wu» Né) cese o + Né) =L (4.9) 
a=1 a=1 
In particular, in the linear example above, we set €; = —1 and č = 1 with n = 2 
(p — 1); then, we obtain the two shape functions using equation (4.7): 
é—1 1 č+1 1 

N = i = == = N: ) pum 1 z 

(O =8@O === =50-9 and ME = HO =le) 
(4.10) 


The global coordinate x can be expressed using the local coordinate č as: 


g = 


2 
N,(£) z; = Ni (ë) af + Na(é) ay = 5 (1 €) a4 5 (1 + €) a. (4.11) 

a=1 

Using equation (4.11), the derivative of the coordinate function is given by: 


=. (4.12) 
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The above derivative between the global and the local coordinates is Jacobian, 
which is represented as je. The derivatives of the shape functions are given by: 
ON, 1 ON, 
=- y @= 1,2, 4.13a 
Je (5 o us 


ON, 10N, 1 ON, 10N 1 
Oz je OE he’ ðr ($06 he 


(4.13b) 


4.1.3 Quadratic Lagrange Element 


For one-dimensional quadratic shape functions, the local coordinates of nodes are 
located at: 


& = —1, Ča = 1, and $3 = 0, (4.14) 


in which the order of the polynomial is p = 3 — 1 = 2 (n = 3). The quadratic 
Lagrange element with three nodes is shown in figure 4.3. 


Element e 
i, Rd att 
3 2 


Node number 1 


Fic. 4.3 — One-dimensional quadratic Lagrange element with three nodes. 


We obtain the three shape functions using the Lagrange interpolation formula 
(4.7) as: 


MO = RG = aR =e 1) 
N(ë) = (5) = eo . Iq 1) (4.15) 
MG) = pO- ERIE 1-2 


If we let the global coordinates for the element be given by zf, z$, and ay (xf and 
x5 are the boundary end nodes), the global coordinate z can be expressed using the 
local coordinate € as: 


x = N,()ay + Ny (£)a5 + .N3(¢) a5. (4.16) 
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The Jacobian is now given by: 


ð 


č 

The value of this Jacobian function is constant only when the coordinate zj for 
node 3 is at the middle point of the element. The derivatives of the shape functions 
are given by: 


- (ccs) (ees) osi petere o2) aam 


js fi (4.18) 


4.2 Two-Dimensional Triangle Element 


4.2.1 Triangle with Three Nodes 


The finite element domain is defined by dividing the two-dimensional continuum 
domain into a set of mesh composed of some triangular elements, as shown in 
figure 4.4a. The linear functions over the triangle with three nodes, as shown in 
figure 4.4b, are constructed by linear polynomials. It is worth noting that this 
two-dimensional triangle element was proposed and used by Courant (1943) and 
Turner et al. (1956) to analyze the plane elasticity problems in the initial develop- 
ment stage of finite element method. 


b, J 


b, 


(a) (b) 


Fic. 4.4 — Two-dimensional domain discretization into triangle elements: (a) triangle mesh, 
(b) triangle element. 
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The approximate linear function on the triangular element with three nodes in 
the Cartesian coordinate system is expressed as: 


qme —[1 xz yl} %2 p =% tort asy, (4.19) 


where $* is the displacement in the element, which can be the displacement u in the 
x-direction or the displacement v in the y-direction; the unknown parameters %1, a, 
and o3 are computed by the displacements of the three vertices of the triangle. 
According to equation (4.19), the displacement on the three vertices can be 
expressed using the corresponding coordinates as follows: 


Qi 1 Ti Yı 01 
95)9—]|l R p a2 P, (4.20) 
$5 1 z yw 03 


where z, and y, (a — 1, 2, 3) are coordinates on the three vertices. The inverse 
matrix of the coefficient matrix in equation (4.20) can be derived as: 


-1 


1 di Uu 1 ay ag a3 
Į d» yo E 2A bi bo bs ; (4.21) 
l B ws € Q G 


where 
Qj = HY3 — BY, b = Y2— YB, C1 = X33 — 15 
= BY — Y3, b=Yy-Y, 0 = 0 — 23 
a3 = Ayo — VY, bs = y1 — yp, C3 = %— X1 


and A = (2b; +b. + 2303)/2 is the area of the triangle. By substituting 
equation (4.21) into (4.20), we obtain 


3 
D ai 
0 lm w ES Qt a4 a a Di Fa 
P aija 2 jjer i eee 
Q9 ?»—|l vv wy 9) EET b b bs > = 9A D b; 
PA l zx y $5 a Q2 B $5 os 
X Coa 
(4.22) 


The geometric interpretation of the parameters b, and c, is shown in figure 4.4b. 

By substituting the expression equation (4.22) of parameters a, into equa- 
tion (4.19), the displacement function in the element expressed by nodal displace- 
ment can be obtained: 
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P = 41 +AT + wy 
3 1 
= a f t ba D a F Ca f a 
2 34 (9s Pat + CaPay) (4.23) 
= 5, (Gat bat cay) Pe 
“= 2A 
According to the above formula, the shape functions can be obtained: 
1 
N.(2, yY) == (Gat baz + cay), a@=1, 2,3. (4.24) 


2A 


The shape function for a = 3 is shown in figure 4.5a, where it can be seen that its 
value at node 3 is 1 and 0 at nodes 1 or 2. More generally, these shape functions 
satisfy. 


1, a=b 
Nalo, Yo) = bab = D P n (4.25) 


where (25, yp) are the local coordinates. Moreover, the completeness condition then 
requires that @°(2y, y») contain any constant c (displacement of rigid body), yielding 


3 3 
(x,y) = So Ns, y)c=c or XO Nela, y)=1. (4.26) 


a=1 a=1 


Using these shape functions, the approximate function of the element is 
expressed as: 


3 
$* - Y Ns y) 66. (4.27) 


a=1 


N, à 


dc 


uu 


Fic. 4.5 — Shape function N; for two-dimensional triangle with three nodes. 
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For example, as shown in figure 4.6a, the two-dimensional plane problem is 
illustrated in which a set of simple linear functions for a parameter @ can be con- 
structed using linear polynomials on a three-node triangle. From a geometric point 
of view, the numerator of the shape function is twice the area of the triangle 
determined by the nodal coordinates Za, y, of the two nodes, and the point (z, y). So, 
in figure 4.6a, we can define 


2A,(z, y) = (aa + bgt + Cay) (4.28) 
and then, the shape functions can be written as follows: 


y, = B9) a= 1,2,3 (4.29) 


This ratio also defines a simple and convenient triangular coordinate system, 
which is called area coordinates (figure 4.6b). We denote these as: 
a=1,2,3. (4.30) 
In addition, it is important to note that these three coordinates are limited by 
the total area such that 
In + Lla+ L3 =1 (4.31) 


The shape function of the three-node triangle can be obtained by using area 
coordinates, as shown in the following formula: 


N,—-L, @=1,2,3 (4.32) 


D 


(a) (b) 


Fic. 4.6 — Triangle element with area coordinates: (a) subareas of triangle, (b) area 
coordinates. 


4.2.2 Higher-Order Triangle Element 


For triangular elements with more than three nodes, the above analysis process for 
three-node triangles also applies. Figure 4.7 shows the first three members of the 
triangular element family. In addition, for the complete set of two-dimensional 
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polynomials, we can use Pascal triangle, as shown in figure 4.8, to represent them. 
Since it is known from the Pascal triangle that the complete order of any polynomial 
can be exactly matched with the number of nodes of the triangle of the same order, 
the shape function of the six-node triangle can be obtained by using the quadratic 
polynomial as follows: 

9 


X2 


Oy (4.33) 


6 


= 0 + TA + YO3 4 ra + Yds d yas 


'This method requires inversing a 6 X 6 matrix to obtain the expression for the 
nodal parameter $$ in terms of the parameters a. The area coordinates used to 
define shape functions can eliminate the necessity for inversing the matrix. 


(a) (b) (c) 


Fic. 4.7 — Triangle element family: (a) linear, (b) quadratic, (c) cubic. 


AAA, 
IPI ITT 


Fic. 4.8 — Pascal triangle. 
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It is simple to directly represent an arbitrary triangle of order M to derive shape 
functions for higher-order triangle elements (Irons et al., 1968). As shown in figure 4.9, 
denoting a typical node aby three numbers J, J, and K corresponding to the position of 
coordinates Lia, Loa, and L3,, then the shape function in terms of three Lagrangian 
interpolations can be written as: 


N, = (L) U (L) Š (L) with I+J+K =M, (4.34) 


where l} , I4, and If are given by equation (4.7), with L4, L2, and L3 taking the place 
of č. 
It is easy to verify the results of the above expressions: 


N,—1 at Lı = Liz, Lə = Lsj, L; Lax, (4.35) 


and 0 at all other nodes. The highest term occurring in the expansion is L/ Lj L4, 
and, as I+ J+ K= M for any point, the polynomial is also of order M. 

Equation (4.34) is valid for quite an arbitrary node distribution in the graph 
shown in figure 4.9, and if the nodal lines are equally spaced (i.e., 1/M), the formula 
can be simplified (Taylor, 1972; Silvester, 1969; Argyris et al., 1968b). The reader 
can easily verify the second-and third-order shape functions as given below and 
deduce the shape functions for higher-order elements. 


(0,0,M) 
2) 


Fic. 4.9 — A general triangular element. 
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4.2.8 Quadratic Triangle Element 
Corner nodes: 
N,—(2L,—1)h, a=1, 2, 3. (4.36) 
Mid-side nodes: 
Ni 4L, Ne=4Algls, Ng =4Ialh. (4.37) 


The concept of the quadratic triangle element was first proposed by 
de Veubeke (1965) and later used in the plane stress analysis problems by Argyris 
(1965). 


4.2.4 Cubic Triangle Element 


Corner nodes: 
1 
N, = 3 3L. —1(3L,—2)L, a=1, 2, 3. (4.38) 
Mid-side nodes: 


9 9 9 
Ny = 5lila(3l1—1), Ns =5Lila(3l2—1), No = 5 Lels(3l2—1) 
(4.39) 
9 9 9 
Ny = 5l -1, Ne=5lsli(3ls—1), No —5laL(3L - 1) 


Internal node: 
Nio = 27 Li Lo L3. (4.40) 


The ‘bubble’ function arises in the last shape function for an internal node with 
no contribution along boundaries. 


4.3 Two-Dimensional Rectangle Element 


4.3.1 Linear Rectangle Element with Four Nodes 


Now we consider the second example of two-dimensional shape functions, the 
rectangle shown in figure 4.10. The rectangular element considered has side lengths 
of 2a in the x-direction and 2b in the y-direction. 

To make it easier to derive the shape functions, we can use a local Cartesian 
system 2’, y' defined by: 


Y —r—a and y =y- y, (4.41) 


44 Basic Theory of Finite Element Method 


Fic. 4.10 — Rectangle element geometry and node numbers. 


where 


1X Le 
% = os z, and y= PD Ya (4.42) 
a=1 a=1 


in which tọ, yo are located at the centre o' of the rectangle, and z, and y, (a = 1, 2,3, 
and 4) are the coordinates of the nodes. For each displacement component, four 
functions are now required to uniquely define the shape function. Its expression is as 
follows: 


qe zo —[1 z y xy] = 0 + T'a + y 0s + a’ y 04, (4.43) 


where $^ is displacement in the element, which can be the displacement u in the zq- 
direction or the displacement v in the y-direction; the unknown parameters o4 to o4 
can be derived using the displacements at each of the four vertices of the rectangle. 
The coefficients a, may be obtained by using equation (4.43) at each vertex node: 


n 1 —-a —b ab 01 
gs {= |1 a —b —ab|]J o 
ef li a 5b ab |] $654) 
ót 1 —a b —ab]| |o 
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By inverting the coefficient matrix and using substitutions, as in equations (4.21) 
and (4.22), we can again solve for a, in terms of the nodal displacements to obtain 


sal Y "ATTEN ES PTR 2T 
e-i- Ge 201) 


4.45 
Iha BV Noth- Nh- Yee e 
' 4 a b P37 4 a b Pa 
which yields the four shape functions: 

/ / / / 

Apei(q2E 1-2), aniio 21-2 

4 b 4 ü b 
(4.46) 


a 
E. Tz V d a! y 
x-i(1 2) 042) m=3(1-2)(1+3) 


To use numerical integration to evaluate the integrals in weak form and obtain 
the rectangular element geometry and the local node numbers shown in figure 4.11, 
we let. 

al $ 
č=— and y= m 
a b 
which are products of one-dimensional Lagrange interpolations using the coordi- 
nates -1 € 6, 4 € 1. Then, the shape functions may be written as: 


(4.47) 


(1-500-9, NM -i0-590-) 


(4.48) 


ALR ORIe 


(1+ (1+), N,—--.0-6(-n 


Fic. 4.11 — Rectangular element geometry and local node numbers: (a) global coordinates, 
(b) local coordinates. 


46 Basic Theory of Finite Element Method 


The shape function for a = 3 is shown in figure 4.12, where it can be seen that its 
value is 1 at node 3 and 0 at nodes 1, 2, and 4. More generally, these shape functions 
satisfy. 


: 1, a=b 
Nile m) = Ôa = Tr um i: (4.49) 


where (č, N,) are the local coordinates. Moreover, the completeness condition then 
requires that @°(€,7) contain any constant c (displacement of rigid body), yielding. 


4 


=X Nén) c=c or Y =1. (4.50) 


a=1 


> 
x 


Fic. 4.12 — Shape function N3 for two-dimensional rectangle with four nodes. 


4.3.2 Higher-Order Rectangle Element 


According to the above illustrations, a systematic and simple way to generate shape 
functions of any order can be obtained by the product of Lagrange polynomials, 
including the two local coordinates (Zienkiewicz et al., 1969; Argyris et al., 1968a; 
Ergatoudis et al., 1968). Specifically, if a node is marked in two dimensions with its 
column number J and row number J, then the following formula can be obtained by 
using the definition of Lagrange interpolation in equation (4.7): 


Na = Ny = YE) Up (4.51) 


in which n and m refer to the number of subdivisions in each direction. For the 
rectangle with four nodes, as shown in figure 4.13a, the mappings of Jand J to node 
number a are shown in table 4.1. 
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TAB. 4.1 — Numbering for rectangle with four nodes. 


Label Node number 
1 2 3 4 
I 1 2 2 1 
J 1 1 2 2 


Figure 4.13 shows a rectangle element with m — n. When m — n — 1, we can 
easily obtain the following result: 


in which ča, Na are the normalised coordinates at node a. 


Fic. 4.13 — Rectangle element family: (a) linear, (b) quadratic. 


4.3.38 Quadratic Rectangle Element 


Table 4.2 shows the mappings of J and J to the node numbers a for quadratic 
rectangle elements, as shown in figure 4.13b. 


TAB. 4.2 — Numbering for quadratic rectangle element. 


Label Node number 
1 2 3 5 
3 3 1 2 2 
J 1 1 3 3 1 2 3 2 2 


The shape functions of quadratic rectangle elements yielding by the products of 
one-dimensional quadratic Lagrange interpolations are as follows: 
Corner nodes: 


1 
Na = 46H (E+ ča) (mm) a@=1, 2, 3, 4 (4.53) 
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Mid-side nodes: 


1 
ča =0, N -5n(1- P) (n+), a=5,7 
(4.54) 


no=0, Ne= (E+E) (1-10), ac 8 


Centre node: 


Ny = (1- 8) (1-1) (4.55) 


4.4 Three-Dimensional Tetrahedron Element 


4.4.1 Linear Tetrahedron Element with Four Nodes 


As shown in figure 4.14, a simple element for three-dimensional problems is a 
tetrahedron with four nodes. 


Fic. 4.14 — Tetrahedron element geometry and node numbers. 


Similarly, a compatible displacement field derived from a complete linear poly- 
nomial expansion is as follows: 


o x —[1 x y z] = 0 + xot + gota + 204, (4.56) 
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where $* is the displacement in the element, which can be the displacement u in 
the z-direction, the displacement v in the y-direction, or the displacement w in the 
zdirection; the unknown parameters a, to o4 may be evaluated in terms 
of the displacements at each vertex of the tetrahedron. The coefficients a, may be 
obtained by using equation (4.56) at each vertex node: 


ĝi laya 9 
$5 _ lam yo 2 02 (4.57) 
$5 1 B Vs % a3 f 
$i l zm Ys A 04 
The inverse is 
-ü 
1 Xi yı ZL [1] a» a3 a4 
lm» y 2 1 bi b bg by 
ae 4.58 
l R Ws B 6V la & cy |’ ( ) 
1 XA YA ZA dy də ds dy 
where 
1 Ty yı zi 
yaga! *? * | (4.59) 
l zs ws % 
l zu Ww au 
with V being the volume of the tetrahedron, and 
D wo 29 l yo g 
àj—det|z ys 2], bb =—det]1 ys 3 
TA Ya 74 lw A ( 46 0) 
l R g l R w l 
C1 = det | 1 B3 «231, dy = —det] 1 % Y3 
1 XA 24 1 TT UA 


with the other constants defined by cyclic interchange of the subscripts in the order 
1, 2, 3, 4. 

The above solution for the parameters a, allows for rewriting the element 
interpolations as nodal parameters: 


Q^ = 03 + tox + oy + 042 


Qaa F box + cay + daz 
(aa ĝ Pay + daĝaz) (4.61) 


(aa + bat + cay + daz) $6 
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This yields the shape functions: 


1 
A 


N.(2; y, z) dat batt Cay + daz), a= 1, 2, 3, 4. (4.62) 

As discussed previously, the tetrahedral element family shown in figure 4.15 is 
similar in properties to those of the triangular family. First, for each stage, we can 
also get complete polynomials in three coordinates. Besides, since the faces are 
divided in the same way as the previous triangles, the same order is obtained for 
the polynomials in two coordinates in the plane of the face and thus ensuring the 
compatibility of elements. In addition, the polynomials contain no redundant terms. 


Fic. 4.15 — Tetrahedron element family: (a) linear, (b) quadratic, (c) cubic. 


'The above form allows the area coordinates to be generalized to the volume 
coordinates represented by L4, Lə, L3, and L4, as shown in figure 4.16, where 


Va 
V ? 
in which Vi, V2, V3, and V4 are the sub-volumes of the tetrahedra 234P, 413P, 421 P, 


and 123P, and V is the volume of the tetrahedron with vertices 1, 2, 3, and 4, 
respectively. These sub-volumes are obtained according to the following formula: 


iS a=1, 2, 3, 4, (4.63) 
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6V, = da + box + Cay + daz. (4.64) 


The shaded region in figure 4.16b is the sub-volume Vj. The volume coordinates 
must satisfy the constraint: 


M- 


ll 
= 


Li= L+ L+ L+ L4 — 1. (4.65) 


T 


'The shape functions of the volumetric coordinate for a four-node tetrahedron are 
as follows: 


N= La, a=1, 2, 3, 4. (4.66) 


{a) 


Fic. 4.16 — Tetrahedron element with volume coordinates: (a) four-node tetrahedron, (b) Vi 
for 234P. 


4.4.2 Higher-Order Tetrahedron Element 


By establishing appropriate Lagrange interpolations similar to equation (4.7), it can 
be seen that the formulae for the shape functions of the higher-order tetrahedron are 
exactly the same as that for triangles. The result is 


N, = (L) U (L2) UE (Ls) (La) with T+ J+K+L=M, (4.67) 


in which M is the order of the tetrahedron. For the following shape functions for the 
quadratic and cubic cases, the readers can verify them directly. 


4.4.8 Quadratic Tetrahedron Element 
Corner nodes: 


N, = La (2La— 1), a=1, 2, 3,4 (4.68) 
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Mid-side nodes: 
N; =4LoI,, Ne =4L03l,, Ny = 4L4Llı 


N c4bA. Dile Diy = beds 09) 
4.4.4 Cubic Tetrahedron Element 
Corner nodes: 
1 
Na = 3 La(3La— 1) (3La— 2), a-1,2,3,4 (4.70) 
Mid-side nodes: 
9 9 9 
N; = z lla -1, N= z liL, -1, M= z Ibis za 
9 9 9 
Na = z PBB — 1), Ng = 5 els(3l1 — 1), Nio = z BlI — 1) 
(4.71) 
9 9 9 
Ny = z Ila (Ls — 1), Ny = z 2243 (3L4 — 1), Ni = z s431 — 1) 
9 9 9 
Nu = z Bla GL; —-1, M; = z Ila GLa -1, Me= z Isla (La — 1) 
Mid-face nodes: 
Ny = 27 Li LəL3, Nis = 2T LAIAL 
17 14243 18 4412 (4.72) 


Nig = 27 L3L4DI4, Nog = 27 Do L3 L4 


4.5 Three-Dimensional Hexahedron Element 


4.9.1 Hexahedron with Eight Nodes 


For a three-dimensional problem, we take the hexahedron shown in figure 4.17 as an 
example for analysis. The derivation of the shape functions for the hexahedron is the 
same as the procedure used for the four-node rectangle. In the present case, it is 
simple and convenient to use a local Cartesian system z', y', z defined by: 


g=2-H, Y=y-y, Z7 =z- %, (4.73) 
where 
1 8 1 8 1 8 
A-3575 W=) Yar = EDA (4.74) 
a=1 a=1 a=1 


in which 25, yo, 2 are located at the centre of the hexahedron, and £a, Ya; Za is the 
coordinates of the nodes. 
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Fic. 4.17 — Hexahedron element geometry and node numbers. 


We may write a polynomial expression for q^ as: 


fy Lay yd da ay Zz] 


^mot-—-|lixyzzxywyzzmxmyz 
o mo [lay y dl MAR (4.75) 


if / / Ld Fa Jd (are Oe i 
HH + T A2 + YO +t zo4-- Vyas +yz as +2 X0 d- X y z etg 


where $* is displacement in the element, which can be the displacement u in the 
z-direction, the displacement v in the y-direction, or the displacement w in 
the z-direction; the unknown parameters o, to ag may be evaluated in terms of the 
displacements at each of the eight vertices of the hexahedron. 

To use numerical integration for the integrals in weak form and obtain 
the geometry and local node numbers of the hexahedron element, as shown in 
figure 4.18, we suppose 


J / z 
, and (——, (4.76) 
a b c 
which are the local coordinates, with the range of — 1 < 7, ¢ < 1. 

The coefficients a, may be obtained analogously as the derivation process of the 
case of the two-dimensional rectangle element. By inverting the coefficient matrix 
and using substitutions, we obtain the eight shape functions: 
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Fic. 4.18 — Hexahedron element geometry and local node numbers: (a) global coordinates, 
(b) local coordinates. 


M-i0-8ü-20-0, &-;ü-00-90-0 
M-i0s00-00-0, N-iüc00-m 040 

; (4.77) 
M= MIU: M=Z0+9 +n) 00 
M-iü0-800400-0, M=Z-9 +n) 049 


4.5.2 Higher-Order Hexahedron Element 


We can describe the equivalent Lagrange-family elements of the three-dimensional 
hexahedron type in a precisely analogous way to that given. Furthermore, we also 
can use a direct product of three Lagrange polynomials to generate the shape 
functions for such elements. By extending the notation of equation (4.51), we now 
have the interpolations: 


N, = Nu = PE Un) BO) (4.78) 
for n, m, and p subdivisions along each side, and with 


Table 4.3 shows the mappings of J, J, and K to the node numbers a for a 
hexahedron with eight nodes. 


Tas. 4.3 — Numbering for hexahedron with eight nodes. 


Label Node number 
a 1 2 3 4 5 6 7 8 
I 1 2 2 1 1 2 2 1 
J 1 1 2 2 1 1 2 2 
K 1 1 1 1 2 2 2 2 
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Figure 4.19 shows a hexahedron element with m = n= p. For m = n= p= 1, 
we obtain the simple result: 


NalE, m O =E EE) Anan) (14+ 600) (4.80) 


in which ča, Na, Ġa are the normalised coordinates at node a. 

The higher-order hexahedral elements were proposed by Zienkiewicz et al. (1969) 
and then elaborated on, particularly by Argyris et al. (1968a). All the comments 
about internal nodes apply here. Figure 4.19 shows the first two elements in the 
three-dimensional Lagrange family. 


Fic. 4.19 — Hexahedron element family: (a) linear, (b) quadratic. 


4.5.8 Quadratic Hexahedron Element 
The shape functions of quadratic hexahedron elements yielding by the products of 
one-dimensional quadratic Lagrange interpolations are as follows: 
Corner nodes: 
1 
Ny = gem (E+E) pem) CHG) a-12,.,8 (481) 
Mid-side nodes: 


1 
Ča = 0, No = qt 1- é) (n 4-17,) (C. a — 9, 10, 11, 12 


u (E+ éa) (1— n?) (C+ ia), a= 13, 14, 15, 16 (4.82) 


(0, Ne té (E+ ču) (N+M) (1-2), a=17, 18, 19, 20 


Ng =0, Na = 
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Mid-face nodes: 


Ča = 7 0, No=5t (=E= n) (Cia), @=21, 22 
Na = a — 0, N,- i (£-&) (1-5) (12 C), a= 23, 24 (4.83) 
Ca = 6,7 0, N,-$n (1-8) (qn) (1-@), a-25, 26 


N-üu-)ü0-9g)(-0O) a@=27 (4.84) 


4.6 Exercises 


4.6.1 One-dimensional element of line with two nodes. 


(1) Using the equations of the shape functions, verify the following: 


diede pM (4.85) 


csse. (4.86) 


(2) Using the equations of the shape functions, compute the following derivatives: 


ON, PN. =| 
oé ’ og , a= l, 


2, (4.87) 


4.6.2 Two-dimensional element of triangle with three nodes. 
(1) Using the equations of the shape functions, verify the following: 


1 =b 
Na(2, y) = Oab m { 0. id b’ (4.88) 


3 


XO Naess (4.89) 


a=1 


(2) Using the equations of the shape functions, compute the following derivatives: 


ON, ON, ON, OPN, ON, 
Ox’? Oy’? Oxy’? Ox?’ Oy" 


a=1, 2, 3. (4.90) 
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4.6.3 Two-dimensional element of rectangle with four nodes 


(1) Using the equations of the shape functions, verify the following: 


Nalo m) = 9a = { 5 ; E $ (4.91) 
4 

>> Nén) = 1. (4.92) 
a=1 


(2) Using the equations of the shape functions, compute the following derivatives: 


ON, ON, ËN, ON, ON, 
oé 3 On a 0én ? og ki On? r 


a= 1, 2, 3, 4. (4.93) 


4.6.4 Three-dimensional element of tetrahedron with four nodes. 
(1) Using the equations of the shape functions, verify the following: 


1 =b 
NS, Yb, 2p) — Oab = ES m "E (4.94) 


4 
5 N,(2, y, 2) = 1. (4.95) 


(2) Using the equations of the shape functions, compute the following derivatives: 


ON, ON, ON, N, N, 
Ox’ Oy? Oz? Oay? Oyz 


4.96 
ON, ON, ON, ON, ON, P NET (590) 
Ozz’ ðr?’ Oy? 82’ Oayz’ DE a 
4.6.5 Three-dimensional element of hexahedron with eight nodes. 
(1) Using the equations of the shape functions, verify the following: 
1, a=b 
Nal Čo tfj, C) = Sab = { 0, a x b (4.97) 
8 
S NEU = 1. (4.98) 


a=1 
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(2) Using the equations of the shape functions, compute the following derivatives: 


ON, ON, ON, ON, N, 
oé’ On’ OC’ DEN’ One 
ON, ON, ON, ON, ON, —— 
aE’ og Op oc OD C 07 


(4.99) 


4.6.6 Summarise the functions of Lagrange interpolation for constructing one-, two-, 
and three-dimensional elements. 


Chapter 5 


Isoparametric Element and Numerical 
Integration 


5.1 Isoparametric Element 


5.1.1 One-Dimensional Isoparametric Lagrange Element 


When the finite element method is applied to one-dimensional problems, elements 
with different shapes and sizes are often used. To unify elements with arbitrary 
length in a one-dimensional problem, a method for mapping elements in the global 
coordinate system to standard elements in the unified local parent coordinate sys- 
tem is developed. Then, all computations in an element are converted to the stan- 
dard element for processing. A representative isoparametric mapping of a 
one-dimensional Lagrange element with two nodes is shown in figure 5.1. 


1 


x I» e————e 
o č 
1 2 1 " 
(a) (b) 


Fic. 5.1 — Isoparametric mapping of one-dimensional Lagrange element with two nodes: 
(a) global coordinates, (b) local parent coordinates. 


The relation between the global and local parent coordinates for a 
one-dimensional Lagrange element is given by: 


a=" NS) ar. (5.1) 


DOI: 10.1051/978-2-7598-2938-5.c005 
© Science Press, EDP Sciences, 2023 
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The use of isoparametric mappings was introduced by Taig (1962). The types of 
parametric mapping that can be considered in the analysis are defined by the 
relation between the coordinate interpolation and the dependent-variable interpo- 
lation. We denote the mapping by: 


a = YO NAE) af = M(E) AHNE USH e Nu IEL (52b) 


a=1 


Once the shape functions are available in terms of the parent coordinates, we 
may immediately use the concept of parametric mapping. We have the following 
three cases: 


(1) Sub-parametric interpolation: The order n of the interpolation for x is lower 
than that order m for u. 

(2) Isoparametric interpolation: The order n of the interpolation for x is the same 
as that order m for u. 

(3) Super-parametric interpolation: The order n of the interpolation for z is higher 
than that order m for u. 


Indeed, we need not always use the same interpolation for coordinates and 
dependent variables; however, an isoparametric mapping is the most convenient 
because nodal coordinates and nodal dependent variables are placed at the same 
physical locations. 

Owing to the transformation of the above coordinate system, the derivatives of 
several functions are involved in finite element computation. This requires switching 
between different coordinate systems. The derivatives of functions on element e with 
respect to the two coordinate systems are related as follows: 


Of Of dx Af CONE), Of. 


dE ^ Or0t dry Ot 7 Os Ps 


1 
here ; — yx”, Ne) vei : ; i 

where j = 55,4 ge is known as the Jacobian for the transformation and by the 

relation defining the local parent coordinates, it can be computed explicitly in terms 

of the local coordinates. 


5.1.2 Two-Dimensional Isoparametric Triangle Element 


The concepts of isoparametric mapping and isoparametric element can be used in 
two-dimensional triangle element elements to obtain a two-dimensional isoparametric 
triangle element. A representative isoparametric mapping of a two-dimensional 
triangle element with three nodes is shown in figure 5.2, in which x and y are the 
global coordinates and L1, L2, L3 are the local parent coordinates. 
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34, OM) 


LR 


AL, LL) AN 
a L = 


(a) (b) 


Fic. 5.2 — Isoparametric mapping of two-dimensional triangle element with three nodes: 
(a) global coordinates, (b) local area coordinates. 


The relation between the global and local parent coordinates for a 
two-dimensional triangle element is given by 


x° = V NI, In, Ls) af. (5.4) 


a=1 


The isoparametric mappings for the coordinate interpolation and the 
dependent-variable interpolation are 


n 3 
z? = M ND L3) zé, O< Li, La, L3 <1 and ML -1 (5.5a) 
a=1 i=1 
n 3 
$° — M NL Ls) @, 0< L, Lo, L3 <1 and ML -1 (5.5b) 
a=1 i=1 


where Li, L2, L3 are the local area coordinates defined above; z^ is the global 
coordinate in the element, which can be the coordinate x or y, and $* is the 
displacement in the element, which can be the displacement u in the z-direction or 
the displacement v in the y-direction. 

By the usual chain rule of partial differentiation, the z and y derivatives of any 
function f may be written in matrix form as: 


Of 3 Of OL; 
Ox i=l OL; Ox 
af 3 Of AL; 
Oy 5L Oy 
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The derivatives of the area coordinates may be obtained by differentiating the 
three equations: 


3 
> Li 
1 i—l 
f = TF = 3 Na (L, Lə, L3) xe . (5.7) 


Using the chain rule, the derivatives of f with respect to z and y may be written 


as: 
OL, OL, 
Ox Oy 
0 0 1 1 1 
of a OL, OL 
ii 0 1 Y| Y) Ys MM 
Oly OL. 
Ox Oy 
where we denote the derivatives by: 
ON. e m "ON, , 
Xj = 3 and Y;— 2. OL, yi. (5.9) 
Solving equation (5.8) yields 
0L, Ol. 
Ox Oy i 
1 1 1| fO 0 B C 
Ol, OL 1 l 1 
ae E = |X X X; 1 0| =a | GO, (5.10) 
P ur Y| Y) Ys 0 1 By O 
Ox Oy 
where 
B= Y- Ya, OC, = X-X 
Baye. Ws (5.11) 
B= Yi- Y, QG=wM-M 


and A= (Xi Bi + X2 Bo + X3 B3)/2. 
'Through the above analysis, the derivatives of the area coordinates derived from 


equation (5.10) are substituted into equation (5.6), and the z and y derivatives of 
any function f can be obtained. 
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5.1.3 Two-Dimensional Isoparametric Rectangle Element 


'The concepts of isoparametric mapping and isoparametric element can be used in 
two-dimensional rectangle elements to obtain a two-dimensional isoparametric 
rectangle element. A representative isoparametric mapping of a two-dimensional 
rectangular element with four nodes is shown in figure 5.3, in which x and y are the 
global coordinates, and č and y are the local parent coordinates. 


Ó 
[52529] 65,55) 


"y 


(a) (b) 
Fic. 5.3 — Isoparametric mapping of two-dimensional rectangle element with four nodes: 


(a) global coordinates, (b) local parent coordinates. 


The relation between the global and local parent coordinates for a 
two-dimensional rectangular element is given by: 


c= 2, N, (é, N) a£. (5.12) 


The isoparametric mappings for the coordinate interpolation and the 
dependent-variable interpolation are given by: 


P=) None, pE5gs1 (5.13a) 
a=1 

$ =X N( né -1s&ns1 (5.13b) 
a=1 


where € and y are the local parent coordinates defined above, z^ is the global 
coordinate in the element, which can be the coordinate x or y, and $* is the 
displacement in the element, which can be the displacement u in z-direction or the 
displacement v in the y-direction. 

By the usual chain rule of partial differentiation, the č and 7 derivatives of any 
function f (£ 7) may be written in matrix form as: 
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Of Ox Oy Of of 
oé Oc O Ó ð 
IM cae er ge (5.14) 
Of Ox Oy Of of 
on On On} V Oy Oy 


The array J is known as the Jacobian matriz of the transformation, and by the 
relation defining the parent coordinates, it can be computed explicitly in terms of 
the local coordinates. In terms of the mapping defining the coordinate transfor- 
mation, we have 


n ON, n ON, 
Pod GEI 
aci 06 gei O6 
J- (5.15) 
)» ON, > ON, 
a, Va ^a. Ya 
a=1 on a=1 on d 


On the other hand, the z and y derivatives of any function f (č, 7) can be 
expressed as: 


of of 
Ox ea oč 
or( J att (5.16) 
dy 0n 


The inverse of the Jacobian matrix is easily obtained from equation (5.4) and is 
given by: 


oy ð 
1| On OE 
J=- , 5.17 
Aloe (5.17) 
ôn OF 
where 
OxOy OxOy 
= det J = —- i 
j= detJ BE On On Be (5.18) 


is the determinant of the Jacobian matrix. 


5.1.4 Three-Dimensional Isoparametric Tetrahedron 
Element 
For a three-dimensional tetrahedron element with four nodes, the isoparametric 


mapping is shown in figure 5.4, in which z, y, and z are the global coordinates and 
Li, Lo, L3, L4 are the local parent coordinates. 
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4 
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3Q 


O 
1 yiz) 


(a) 


Fic. 5.4 — Isoparametric mapping of three-dimensional tetrahedron element with four nodes: 
(a) global coordinates, (b) local volume coordinates. 


The relation between the global and local parent coordinates for a 
three-dimensional tetrahedral element is given by: 


x° = V NI, Lo, Ds, L4) £$. (5.19) 
a=1 


The isoparametric mappings for the coordinate interpolation and the dependent 
variable interpolation are given by: 


a 


a: 


4 
Na(L1, Lə, Ls, L4) z, O< Li, La, Ls, L4 <1 and ŅX`Li=1, (5.20a) 


n 
=1 l 


n 4 
@° = M ND Ls, L4) @§, O< La, In, Es, I4 <1 and ML —1, (5.20b) 


a=1 i=1 


where L1, L2, L3, L4 are the local volume coordinates defined above, z^ is the global 
coordinate in the element, which can be the coordinate z, y, or z, and @° is the 
displacement in the element, which can be the displacement u in the z-direction, the 
displacement v in the ydirection, or the displacement w in the z-direction. 

By the usual chain rule of partial differentiation, the x, y, and z derivatives of any 
function f may be written in matrix form as: 


of 1 Of OL; 
Or i=l OL; Ox 
of 1 Of OL; 
4 y= 21 
Oy 2L, Oy pos 
Of 4 Of OL; 
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The derivatives of the volume coordinates may now be obtained by differenti- 


ating the following four equations: 


pi 


RC 8 
^ d ^ 


>> Na(l, Lo, Lg, L4) z 


a=1 
n 


D NL, Lə, La, L4) y 


a=1 
> NL, Lo, La, L4) Zz, 


(5.22) 


Using the chain rule, the derivatives of f with respect to z, y, and z may now be 


written as: 
0.0 0 
Of of of 100| _ 
Ox Oy Oz} |0 10| - 
0 0 1 


where we denote the parent derivatives by: 


X,= aE Y; = ae 


Solving equation (5.23) yields 


ðr Oy Oz 
Oly Oly Oh| [i 4 
Or Oy Oz X; X 
Ola ð Lz ð Lz i Yi Y 
Or Oy Oz 4 2 
ðr Oy Oz 
where 
1 Yo & 
B, = — det 1 Y3 Z3 ; Ci 
1 Yi Z 


1 1 1 1 Ol OLy Oly 
X; X; X; Xs! | Ox Oy Az 
Y Y Ya Yı ð Lz ð Ls ð Ls 
ZA 2 B ZA Oy Oz 
OL, Ol, OL, 
Or Oy Oz 
ONa ye ye "ON, , 
and Zi-— OL, Za 
i] “fo 0 0 B 
X% X 100 i 
Ys Yi 010| 6V] Bs 
Z Z 00 4 B, 
1X» 2 1 
= det|1 Xs Z3 ; D; = — det 1 
1 X4 ZA 1 


C 
C5 
C3 
C4 


X» 


X4 


(5.23) 


(5.24) 
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with the other constants defined by a cyclic interchange of the subscripts in the order 
1, 2, 3, 4, and 


1 X Y A 
ule lt Be Ex 
1 X4 Yn A 


Through the above analysis, the derivatives of the area coordinates derived from 
equation (5.25) are substituted into equation (5.21), and the z, y, and z derivatives 
of any function f can be obtained. 


5.1.5 Three-Dimensional Isoparametric Hexahedron 
Element 


For a three-dimensional isoparametric hexahedron element with eight nodes, the 
isoparametric mapping is shown in figure 5.5, in which z, y, and z are the global 
coordinates, and č, y, and ¢ is the local parent coordinates. 


y 


(a) (b) 


Fic. 5.5 — Isoparametric mapping of three-dimensional hexahedron element with eight nodes: 
(a) global coordinates, (b) local parent coordinates. 


The relation between the global and local parent coordinates for a 
three-dimensional hexahedral element is given by: 


= s Nm) at. (5.28) 
a=1 


The isoparametric mappings for the coordinate interpolation and the 
dependent-variable interpolation are given by: 


r= Nén O i -1<é,9,6<,1 (5.29a) 
a=1 
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ĝ = N,(é,n, C) $5; —1<č,ņn, <1 (5.29b) 


a=1 


where €, 7, and ¢ are the local parent coordinates defined above, z^ is the global 
coordinate in the element, which can be the coordinate z, y, or z, and @° is the 
displacement in the element, which can be the displacement u in the z-direction, the 
displacement v in the ydirection, or the displacement w in the zdirection. 

By the usual chain rule of partial differentiation, the £, y, and ¢ derivatives of any 
function f (£ 7, €) may be written in matrix form as: 


of Ox Oy Oz] (ðf Of 
aé OE 0€ OE] | Ox Ox 
Of Ox Oy Oz Of Of 

= = i 5.30 
ðn (| Oy Oy 9| J p (| dy Ven 
of Ox Oz Oz| | Of of 
a at a at] \ az Oz 


The array J is known as the Jacobian matriz of the transformation, and by the 
relation defining the parent coordinates, it can be computed explicitly in terms of 
the local coordinates. In terms of the mapping defining the coordinate transfor- 
mation, we have 


"ON, LƏN ON, 
Ta a z Ža 

Lan Do” Mw 

n ƏN LƏN LÔN, 

J = Ta a Za 5.31 
» On m on ^ M on Oe) 
"ON, LƏN LƏN, 

Ta m a m Za 
ba^ ha^" ba 


On the other hand, the z, y, and z derivatives of any function f (é, 7, &) can be 
expressed as: 


of of 

Ox OE 

Ofl _ 7-1) 9f 

7 =J n[ (5.32) 
of of 

Oz at 


The inverse of the Jacobian matrix can be easily obtained from equation (5.31), 
and is given by: 
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Jt = |b b b (5.33) 
J [4] C2 C3 
where 
n ON, n ON, n ON, n ON, 
UA a ^a “a a he a fa 
2 On i 2 On m On m on 
a — 09 V aw. nan, |^ 47799 , an n ON 
“Ya 7 da T ds = Ba 
Sa” dar? La” LX 
(5.34) 
n ON, əN, 
M p > ðq ^ 
a RES n ON, 
ae Hor” 


with the other constants defined by the cyclic interchange of the subscripts in the 
order 1, 2, 3; 


j= det J (5.35) 


is the determinant of the Jacobian matrix. 


5.1.6 Requirements of Isoparametric Element 


As noted above, the mapping from parent to global coordinates in an element must 
always be single-valued. That is, each point in the parent coordinates should have 
unique global coordinates. The parametric mapping is controlled by the placement 
of nodal coordinates for each element. The condition that ensures a one-to-one 
mapping is that the Jacobian determinant j (computed by the chain rule) is 
positive: 


j20 (5.36) 


For lower-order elements, one can establish criteria to ensure that equa- 
tion (5.36) is satisfied everywhere in the element. However, as the order increases, 
this becomes more difficult. To ensure that the condition is satisfied, the following 
must be checked: 


(1) Placement of edge nodes to make j positive. 
(2) Placement of face nodes to make j positive. 

(3) Placement of internal nodes to make j positive. 
(4) All vertex angles are less than 180°. 


Checking that the value of j is positive at each node of an element is sufficient to 
ensure that the above conditions are satisfied, and thus the specified parametric 
mapping produces a valid element. 
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5.2 Numerical Integration 


As the shape functions and derivatives are known, the next step requires the eval- 
uation of integrals in the weak form to obtain the element matrices. 


5.2.1 One-Dimensional Integration for Lagrange Element 


Because the shape functions are expressed in terms of the local parent coordinate, it 
is also convenient to express the integrals in terms of the local parent coordinate. 
The differential is given in terms of the Jacobian as follows: 


Ox 


de ae 


dé = j(£)a€. (5.37) 
Thus, the integrals may be expressed as: 


f toa f focos (5.38) 


where f(¢) = f(z(£)). 


If the integrand is complex, it is difficult to compute the above definite integral. 
In another manner, numerical integration is used to approximate the integrals. An 
efficient formula for polynomial forms is the Gauss quadrature or Gauss-Legendre 
quadrature, which is defined on the interval —1 < € € 1. The integral is replaced by 
the following form: 


[ Hae Dieu, (5.39) 


with an appropriate choice of integration points €; and weights wj. The Gauss- 
Legendre quadrature formula is exact in the case of polynomials of degree at most 
2n-1. Table 5.1 provides typical locations of the integration points and weights. 


TAB. 5.1 — Gaussian integration points and weights for 
one-dimensional Lagrange element. 
n ej Wj 
1 0 2.0 
2 +1 /V3 1.0 
3 +v0.6 5/9 
0.0 8/9 
+V34+V48 1 = 1 
7 2 3/4.8 
i +V3- v4.8 : + : 
T 2. 3/48 
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Example 5.1. Stiffness for a one-dimensional element with two nodes by numerical 
integration. 


The shape functions and Jacobian for a one-dimensional element with two nodes 
are 


Ox he 


1 1 l 
M-30-0, M=5(1+9, j= (5.40) 
The stiffness for a two-node element can be computed as follows: 
; ON, i ON, 0C 
"YBa Ap SON: BN BE Ox | p [ONE ƏN 0E]. 
K? = Ox E, 1 elav >i x n ou ON2 OS d 
| ae E: Ox a) amo ("| BE ae ros] 05 
Ox 0€ Oz 
ON, 1 
ll are zx 
= ot j(&) E E 1 ON, 1 IL j 
Jan P (ae xe “ae xe 
06 jC) 
ON, 2 
: h. 1 
_ OE he ON, 2 ON, 2} 1 = Ef 1 -1 J g 
a ƏN 2 z| ðE he 8€ he g hed = ap E 1 25 dé 
OC he 


which is identical to that obtained previously. 

Numerical integration is performed using the one-point quadrature formula 
because the shape-function derivatives are constant. The numerical integration is as 
follows: 


1 1 
i. dé =X f(&)w; = f(O)m 21x 2-2, (5.41) 
B mq 


which is obviously the correct answer. For this problem, quadrature was not 
required; however, for higher-order elements, quadrature greatly simplifies the 
calculations. 


Example 5.2. Stiffness for a one-dimensional element with three nodes by numerical 
integration. 


The shape functions and Jacobian for a one-dimensional element with three 
nodes are 


M-i-1, M= E+), M=1-2, (5.42) 


30-5 (6-5) ah (tg) ne opere rm. Gam 
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For the case where the coordinate for node 3 is centred between nodes 1 and 2, 
the derivatives are linear in €, and the Jacobian is constant and equal to Ih. 
'The stiffness for a three-node element involves the integral: 


oM 
Ox 


m T ON» z| ON» 8 | a 
0 Ox ðr Ox Oz 
ON; 


Ox 


ON, OC E p OG ƏM, oé ON30E 


OC Ox oč Ox OE se iat 


BER OF TO c sg PO’ 


OE jel) 


1 
| 5 hedé 


aN, 2 p[OM2 952 9N2 
"| 06 he OER, OE he 


ON, 1 
96 j«(&) 
i ON. 1 B| 1 ON. 1 ON, 1 


ag, f (e-a) ( ;) Ga (: ; 20) 
vil (e-3)(e+5) (3) (ese 
(: -s)cxo (e 5) (—2€) (-28y? 


In this case, the exact stiffness is recovered using a two-point quadrature formula. 
For example, Ki, is 


E (C3 0+ (4-4) >) (544) 
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It can be seen that this result is exactly the same as the result of direct inte- 
gration, as shown below: 


2E, f! J^. 2E 
ku =" f (6-3) a= EL GE 


3 2 x 
2E (C 6 + (5.45) 
he \3 2 4/|4, 
TE, 
— 3h. 
'The reader can verify that the final result is 
E, 7 1 -8 
K*-— 3h 1 7 —-8}. (5.46) 
^[—-8 —8 16 


If node 3 is not centred between nodes 1 and 2, the Jacobian is not equal to ih. 
and should be regarded as a part of the integrand. 


5.2.2 Two-Dimensional Integration for Triangle Element 


The numerical integration formula in terms of the area coordinates is given by 
(Anderson et al., 1968; Hammer et al., 1956; Radau, 1880): 


1- Lı 
if n fa, L2, L3) dad Ly «Mr Lij, Loj, Laj) wy, 


EEEE 


(5.47) 


A series of necessary integration points and weights are listed in table 5.2 (Lyness 
and Jespersen, 1975; Cowper, 1973; Abramowitz and Stegun, 1965). 


TAB. 5.2 — Gaussian integration points and weights for a two-dimensional triangle element. 


Order Figure Integration points Area coordinates Weights 
Linear 1 1 
(n=1) : 2g ; 
11 0 1 
a 9' 9? 9 
Quadratic a : 
" 1 0 1 1 
(n = 3) 2 $ V3 2 3 
r 0 1 1 1 
2" 2 3 
211 1 
° 8 B B 3 
Quadratic 
b 121 1 
(n — 3) LN. pg 3 
" 112 1 
6'6'3 3 
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Tas. 5.2 — (continued). 


Order Figure Integration points Area coordinates Weights 
n 111 27 
3'3'3 48 
. 25 
Cubic LN b 0.6, 0.2, 0.2 E 
(n=4) 2 
IIN c 0.2, 0.6, 0.2 z 
d 0.2, 0.2, 0.6 29 
48 
a a L l 0.2250000000 
3° 3°3 
b on, Bis By 0.1323941527 
a c Bi, o2, fi 0.1323941527 
t 
Pm : A d Bis Ba, % 0.1323941527 
= EAA e 99, Bo, Bo 0.1259391805 
f Bo, o2, Bo 0.1259391805 
g Bo, B», %2 0.1259391805 
with 


ot; = 0.0597158717 
B, = 0.4701420641 
a2 = 0.07974269853 
By = 0.1012865073 


5.2.8 Two-Dimensional Integration for Rectangle Element 


Even though integrals for rectangular four-node elements may be calculated quite 
easily, the integration becomes considerably more difficult for axisymmetric and 
mapped elements using the parent coordinates. Thus, in general, it is best to use 
quadrature (numerical integration) as a general computational tool. To use isopara- 
metric mappings, it is first necessary to change the domain over which the elements are 
integrated. This is a standard procedure that is taught in integral calculus. Thus, for a 
two-dimensional problem involving rectangular elements, the integral over the ele- 
ment domain is converted to an integral over the parent domain —1 € ¢, y € 1. 

The infinitesimal area elements in the global and local coordinate systems are 
related as follows: 


dady = dé x dy 
ya Op ic y, , 
= (apti aga) x (Gantt yn) 
f 05 4 OY Oy 
= (Beit Sel) IC (5.48) 
Ox Oy 
8E BE 


oé 
= dédy = jdčdn 
Ox Oy 


an Om 
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Using the isoparametric form, the integral is transformed as: 


[ L F(a, y)dady — T [senate nac, (5.49) 


where j(€,7) is the determinant of the Jacobian matrix described above. With the 
above choice for the domain of the parent element, it is now convenient to evaluate 
integrals using Gaussian quadrature. Accordingly, we can integrate into two 
directions using 


N M 
A L FE nj (E med = V^ YO Fém de Ems My) Wm tn, (5.50) 


n=1 m=1 


where m denotes the quadrature points in the c direction, and n the quadrature 
points in the 7 direction. The location and weight of the quadrature points in each 
direction are given in table 3.1 for one-dimensional integrations. 


Example 5.3. For a two-dimensional rectangle element with four nodes, compute 
numerically the integral fi fi (£ +n? )dédn. 


In this case, the integral is recovered using a two-point quadrature formula in two 
directions: 


1 1 2 2 
l I (E F n? )dédn x > 5 Ems Nn) Wm Wn 
—14-1 


n=1 m=1 
= g(&i,m)wiwi + 9(€1, N2) wi w2 
+ g(č2, N1) wow; + g(&o, N2) uou» 


=o(-.-5) x1x Les) x1x1 
(8-2) xix t ) x1x1 


(5.51) 


C2| oo 


It can be seen that this result is exactly the same as that of direct integration: 


1 pl 1 pl 1 pl à 4 8 
/ / (€ -ip)dédn = J J &dédg + J / n'dédy= + = —. (5.52) 
-14-1 -14-1 -1J-1 3 3 3 


5.2.4 Three-Dimensional Integration for Tetrahedron 
Element 


The numerical integration formula in terms of the volume coordinates is given by: 
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d pe Ae n 
I x D J f, Le, L3, L4) dLsdLod La ~ 3 (By os, Laj, Laj) wj, 
o Jo 0 


j=l 
Sia. i 
(5.53) 


Furthermore, table 5.3 presents some integration points and weights. 


TAB. 5.3 — Gaussian integration points and weights for three-dimensional tetrahedron 
element. 


Order Figure Integration points Volume coordinates Weights 

Linear 11 1 1 
a 4 Arar 4 1 

(n= 1) 4'4'4'4 
1 
a X, B, B, B 4 
h B. o, B, P i 
> %, P, Vi 
Quadratic i 4 
1 
(n — 4) c b, B, «, B Fi 
1 
d B, B, B, X 4 

æ = 0.58541020 
p = 0.13819660 

1111 4 
a Ak At AD A m" 
4'4'4'4 5 
i 1111 9 
, 2'6'6'6 20 
Cubic 1111 9 
(n — 5) MEA : 62866 20 
d 1 1 1 1 9 
6 6 2'6 20 
1111 9 
: 66602 20 


5.2.5 Three-Dimensional Integration for Hexahedron 
Element 


Integration for hexahedron elements is performed as in the case of quadrilateral 
shapes. For a three-dimensional problem, the quadrature over the element domain is 
converted to an integral over the parent domain —1 < é,7, < 1. 
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The infinitesimal volume elements in the global and local coordinate systems are 
related as follows: 


dadydz = dé - (dy x dC) 


Oy 
-(& "dei E apu). 
Oz, ., Oy, , 9? Or... Oy 
| (ean | ay td Zane) x (Fac | ae sack) | 
= (Fi OV ae Jat 


oí oc oë 
Ox, Oy, Oz Ox, Oy, Oz (5.54) 
SC E E WoT EE 
Ie Tanta je * (S^ at! * at Jae 
Ox Oy Oz 
oč O€ OE 
| |Oxr Oy Oz Lais us 
=| ðn By = dd 
ðr Oy Oz 
ot OF oç 


Using the parametric form, the integral is transformed as follows: 


1 


where j-(¢, n, C) is the determinant of the Jacobian matrix for the three-dimensional 
hexahedron element. The numerical integration formula is 


J l i l / Men 


N M 
ddA Flëm Nn EDI (Ems Nn C1) Wm Wn Wi 


n-lm 


(5.56) 


M- 


ll 
un 


where m, n, and / denote quadrature points in the £, 7, and ¢ directions, respectively. 
The location and weight of the quadrature points in each direction are given in 
table 3.1 for one-dimensional integrations. 


Example 5.4. For a three-dimensional hexahedron element with eight nodes, com- 
pute numerically the integral fon ee Jj (E +n? + C)aédgac. 


In this case, the integral is recovered using a two-point quadrature formula in 
three directions: 
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2 2 2 
Pf fe (E i? +0) dédnde Qo »X (Enns ins C1) Win ts uy 
El #=1 m=i 


= g(& m, 01) wr w w + (61, mi, C2) wi w w 
+ 9(S1, M2, C1) wi wo wy + 9(C1, M2, C2) wi wo We 
+ g(&2, 1, G1) wy wi wy + g(&o, N1, C2) uo Wy We 
+ g(a, N2, G1) uy wo wy + g(&o, No, C2) wo wo w» 
1 1 1 1 1 
"(Cou ME M V8 Và 


1 1 1 1 1 
; ; x1x1x1l y , xixixi 
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C V3 z) Ge We 5x) 
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=8 


It can be seen that this result is exactly the same as that of direct integration: 


LL ? +n? + C)dédgdt 
= L 2 J  P'acanit + J l / J l nedédndt + / i / i Caganat (5-58) 


oO ud 
3 3 3. 


5.2.6 Required Order of Numerical Integration 


For one-dimensional problems, with the appropriate choice of points ¢; and weights 
wj, the numerical integration formula is exact for polynomials of degree at most 
2n-1l. For two-dimensional and three-dimensional problems, the number of inte- 
gration points in each dimension need not be the same, but for convenience, the 
same number is often selected. 


5.9 Exercises 


5.3.1 Based on the parametric mapping defined by the relation between the coor- 
dinate interpolation and the dependent-variable interpolation for a one-dimensional 
Lagrange element, describe sub-parametric, isoparametric, and super-parametric 
interpolation. 

5.3.2 Consider the two-dimensional rectangle isoparametric element as shown in 
figure 5.6. 
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(a) Write the expression for an isoparametric coordinate mapping in this element. 

(b) Determine the location of the local coordinates c and y that define the centroid 
of this element. 

(c) Compute the expression for the Jacobian matrix J for this element and evaluate 
the Jacobian matrix at the centroid. 

(d) Compute the derivatives of the shape function N3 at the centroid. 


Fic. 5.6 — Two-dimensional rectangle element. 


5.3.3 Compute the stiffness K;3 for a one-dimensional element with three nodes as 
follows 


(a) By numerical integration using one, two, and three integration points. 
(b) By direct integration. 


5.3.4 For a two-dimensional rectangle element with four nodes, compute numerically 
and by direct integration the integral: 


i l T l (ë+ n)dédn. (5.59) 


Chapter 6 


Finite Element Computation Scheme 
of Elasticity Problems 


In this chapter, we consider two- and three-dimensional elasticity problems 
(Timoshenko and Goodier, 1970) and describe all the necessary steps for the finite 
element analysis of general solids, including the weak form, finite element method 
for solving elasticity problems, global assembly using the element location vector, 
and treatments on boundary conditions. 


6.1 Weak Form for General Elasticity Problems 


We begin by considering the weak form of three-dimensional problems in elasticity. 
The weak form of the equilibrium equations for linear elasticity may be written as: 


G(w,u, o) = [we +S'a)dQ = 0. (6.1) 


Expanding the equations for the three-dimensional problem in Cartesian coor- 
dinates yields 


Oo, + Otys L Ota 


E b. Ox Oy Oz 

G(w,u, 0) = i Wy by $+ a $ x + ors dQ=0. (6.2) 
Ww b. Ot, F Oty, OG, 
Ou Oy Oz 


For the three-dimensional problem in Cartesian coordinates, the infinitesimal 
volume element is given by dQ =dadydz. Integrating the terms involving the 
derivatives of stresses and using Green’s theorem yields 
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where dI = dA is the infinitesimal area on the boundary surface. To illustrate the 
transformation process in the above formula, the derivation and transformation 
process for the first term is presented in detail: 


ð P 2b Yb Tp ð ? 
| w, aQ = f / i w < dadydz 
o  Ór za Jum Ja ôx 
2b Yb Th o Zb Yb 
2 -f J ue o,dadydz + / T [WuFx] |r) — [WuPx]| ;- ad vd z 
Za Ya Ta Ox Za Ya 
2b Yb Th Owu Zb Yb 
mE ji I j - c, dzdydz + f / [WuF2Ne]|e=b + [wo sna] |; dE 
Za Ya Ta Ox Za Va 


o u "t " 
= — H o,dQ + f f Wy Fn 
o Ox Za a 


m t = o,dQ+ | w,c,n;dF (6.4) 


A key transformation in the above formula is 


2b Yb 2b Yb 
/ / [Wu s] |; — n = [Wu x||e—-adydz = J J [Wu 4] z= + [WuFeNe||x—adT . (6.5) 
Za Va. Za Ya 
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This equation can be proved using the projection of the boundary surface of the 
three-dimensional solid on the boundaries + = b and x = a in the coordinate plane 
y-o-z, as shown in figure 6.1. The following geometric relation holds true for each 
infinitesimal boundary surface: 


dydz=n,dI, onthe boundary z = b (6.6a) 


—dydz=n,dI', ontheboundary x = a (6.6b) 


dydz=cos f dI 
=n. 


dydz-cos(n-a)dI" 
=cos a d —n,dI" 
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FiG. 6.1 — Geometric relationship for the integration on the boundary surface of the 


three-dimensional solid. 


Subsequently, we note that 


f) 
de i : 
0 By 0 
0 0 2 Wu 
a ə d Wy p = Sw = £y. (6.7) 
Oy às DA 
ð oð 
0 ð 
sc 5 
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Thus, this is merely the strain variation expressed in terms of displacement 
derivatives. 
We also note that the boundary integral term may be written as: 


Ox 
Wy T n; 0 0 m, 0 m Ei 
Í Wy 0 n, O m mn O0 a dr 
T| wy 0 0 m O m, n g (6.8) 
Cig 
Uu 


= fi w'G'odI = fy wtar 
r 


where t = G'e defines the boundary traction vector. If we split the boundary into 
two parts in which traction and displacement boundary conditions are defined, we 
may insert the specified traction boundary conditions as: 


» wltdr = J w'tdr+ | w'tdr = li wtdF4 | wtar. (6.9) 
Tr D, T; Ty Tr; 


Generally, we impose the displacement boundary condition u = u directly and 
set w= 0 on [,. Thus, the weak form for the linear elastic equilibrium equations 
may be written compactly in matrix form as: 


G(w,u, 0) = 1 [w'b —sLe]dO - | w' tdr — 0. (6.10) 
Q T, 


The weak form of linear elasticity may be written completely in terms of dis- 
placements to compute strains using geometric equations: 


G(w,u) — n w'bdQ— | &LDsdO-- | w'tdr =0, (6.11) 
Q Q T, 
where it is also understood that the strains are expressed in terms of displacement 
derivatives. This is called the displacement method or irreducible form of the problem, 
as it is only necessary to define the distribution of the displacement field to solve a 
problem. 


6.2 Finite Element Method for Solving Elasticity 
Problems 


This basic structure is fundamental to all finite element methods. As the strains 
involve the first derivatives of displacement components, the approximation to these 
components needs only be continuous in Q. A similar result holds for the relation 
between the strain variation and the displacement variation. Thus, in all multidi- 
mensional elasticity problems, we only require Cp continuous functions of the type 
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presented in the previous chapter to approximate the displacement and virtual 
displacement. Accordingly, we assume that the displacements are given in each 
element by: 


wr w=” N(x) dr, (6.12a) 
a=1 

u x ù= Ni(x) à; (6.12b) 
b=1 


where N,(x) and N;(x) are the shape functions, and wt and w, are the displacements 
at the nodes. 
Inserting the E EDDA for displacements into the weak form yields 


G(w,u) ~ G(w, à) = 2 G°(w, à) = (6.13) 


The expression for each element domain Q, is given by: 


G°(w, à) = | i bdQ— | eL DedQ+ I w' tdT, (6.14) 
Q, Tet 


A 


in which the strains are expressed as: 


x Ti (6.15) 
&— M S(N,) i; =X Biù; 
b=1 b=1 
where B, is the three-dimensional strain-displacement matrix: 
ON, 
Ox as 9 
0 l o 
dy 
ON, 
0 0 E, 
B, = ON, ON, 7 . (6.16) 
- 0 
Oy Ox 
0 ON, ON, 
Oz Oy 
ON, 0 ON, 
Oz Ox 


Evaluating the weak form for an individual element yields 


- »» »» wi Ks E» weft = 2 wT (£c Kafe ] (6.17) 
a=1 a=1 b= 


a=1 b=1 
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where the element stiffness and load matrices are computed as: 


k= | B'DB,dQ, fe = n N,bdQ+ I N,tdT. (6.18) 
eA Qe Tet 


In this form, we again obtain a standard discrete system in which the total 
problem is assembled from the contribution of each element. After the assembly of 
all elements, the weak form for the approximation is 


G-Y wl (>: Kaĉo 4) =0 (6.19) 
a=1 b=1 
with 


Kay = S °K‘, fa = If (6.20) 


and as w, is arbitrary, the governing equation of the problem becomes 


X Kavtts =f, forall a (6.21) 
b=1 


We can derive the stiffness equations: 
Ku=f, (6.22) 


which, after the boundary conditions are formally inserted, yields the solution for 
nodal displacements: 


à —K f, (6.23) 


in which K^! is the matrix inverse of K. 
Once the displacements are known, the approximation for strain and stresses in 
the element may be computed as: 


n n n n 


&$— V s(N)ü,— M Bå; ê= M DS(N,) a; — M DB ùi (6.24) 


b=1 b=1 b=1 b=1 


6.3 Global Assembly from High-Dimensional Elements 


The global assembly of two- or three-dimensional elements is exactly the same as 
that in the previous one-dimensional problem according to the element location 
vector. For expository purposes, a typical two-dimensional example is provided to 
illustrate the global assembly process and computation procedure. Some steps in the 
computation are presented in detail; experienced advanced readers may skip them. 
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Example 6.1. Solutions of finite element method for two-dimensional plane strain 
problems. 


A two-dimensional plane strain model is shown in figure 6.2. The dimensions of 
the model are a= 30 in the horizontal x-direction and 6 = 10 in the vertical 
y-direction. The material parameters are Young's modulus E of 1 x 10° and 
Poisson's ratio v of 0.3. The boundary conditions are that the upper left and upper 
right vertices, as well as the lower parts, are fixed. The vertical concentrated load 
F acting on the middle domain is —1 X 106. 


Element number: (D and © Node number: 1-6 


Fic. 6.2 — Element and node numbers of two-dimensional problem. 


Then, for each element in the local parent coordinate system, a standard 
two-dimensional rectangle element with four nodes is used, as shown in 
figure 6.3. 

The relation of the local nodes to the global node number is indicated in 
table 6.1. 

According to the node numbers of each element, the corresponding element 
location vector can be obtained to determine the relation between local and global 
node locations as follows: 


à-(213 4H, P={4 3 5 6M. (6.25) 


'To obtain the global stiffness matrix, the stiffness matrix of each element should 
be computed according to the formula introduced in the previous section: 


K= | Bi DB,dQ. (6.26) 


e 


By substituting the corresponding parameters of this model, the stiffness 
matrices of the two elements shown below can be obtained. Here, the final results are 
given directly: 


K°= 


Ki, 
Ka 
Kj 
Ki 


Ki; 
K5, 
K5» 
Ki 


Ki 
K3; 
Kj; 
Ki 


Kj, 
K5, 
Kj, 
Ki 


491480.77 240384.60 
158546.95 


Symmetric 


—42735.03 
48076.92 
491480.77 


—48076.92 
—630341.82 
—240384.60 

158546.95 


—245726.48 
—240384.60 
202991.45 
48076.92 
491480.77 


—240384.60 

—379273.48 
—48076.92 
251068.35 
240384.60 
758546.95 


—202991.45 
—48076.92 
—245726.48 
240384.60 
—42735.03 
48076.92 
491480.77 


48076.92 
251068.35 
240384.60 

—379273.48 
—48076.92 

—630341.82 

—240384.60 
758546.95 


(6.27) 
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fg 


Fic. 6.3 — Two-dimensional rectangle element with four nodes in local parent coordinate 


system. 


TAB. 6.1 — Local-to-global node numbering for two-dimensional problem. 


Element number 1 2 

1 2 4 

Local node number 2 1 3 
3 3 5 

4 4 6 


'To explain the derivation of the above element stiffness matrix, the computation 
of the first coefficient 491480.77 in K 3 for element 1 is presented in detail. First, the 
matrix forms of the element stiffness matrix are expanded: 


Ka = f B;DB,dO = f ST (N,)DS'( N;)dQ 
Q: [0 


e 


ON, ON, 
Ox pe on Oy 
ON, ON, 
) Oy j Ox 
am o 
=f (1 =o v) U v 0 Ox dQ 
E v (1— v) v 0 0 a 
d v — (1-v) 0 0 a (6.28) 
0 0 0 (1 — 2v)/2 ON, ON, 
Oy Ox 
ONAN yy ON, ƏN, 
Ox Ox Ox Oy 
ON, ON; ON, ON, 
—— (1 — 2v)/2 t (1 — 2v) /2 
n Fal 29/ a acr 
o, ) d ƏN ðN, ON ON. j+ 
Oy Or Oy dy i 
ƏN, ðN, ƏN, ðN, 
ds a 4-20 F a4 20)/2 
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where d = (1 + v)(1 — 2v). 

It can be seen that the above formula involves the derivatives of the shape 
functions with respect to the global coordinates z and y, which can be converted to 
the local coordinate system č and y as follows: 
ON, ON, Oy Oy] ( ƏN, 

_ pi} S04 i| Op a|] a 
ON, ON, j | Óx Ox | ON, 
Oy on On OC on 


n ON, n ON, n ON, n ON, 
( ( Wai EZ ts) ( 324 En 2 = ( Dong an ws) ( Se Ot «) 
n ON, n 
e On Ya Da oé Ya 


n ONa n 
=a >a Ya =1 ae Ta 
2523 On m = 


n ON, n ON, n ON, n ON, 
| ( M E ws) ( des EN 2 e ( DAMES EN 2 ( VAR Ed 2) 
n ON, ON, n ON, ON, 
( vem on 2 ðE ( bacs OE 2 On 


n ON, ON, » ON, ON, 
-(zz ôn ws) aé «(3 Er ws) an 


(6.29) 


where 
oF | OxOy OxOy | n ON, n ON, n ON, 
Tod ean Onde” (es) (3X me) - (Yes) 


n ON, 
ee oé in): 


Here, the following shape functions are used: 


1 1 
NM-40-9: 0-9, NM-40-9- (+n) 
6.30) 
1 1 ( 
N — qr On. N,-q4üc5-0-m 
'The derivatives of the shape functions are 
ON 1 ON... 1 ON; 1 ON, 1 E 
ðN; dou ON; lj. ON, 1 ON, ls E 
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Moreover, the global coordinates of the nodes are 
(n; y) — (0, —10), (x, y») = (0,0), (23, y3) = (15,0), (24, yi) = (15, 2-10). (6.32) 


Using the derivatives of the shape functions and these coordinates, we obtain 


S 
II 


"AN, 1 1 1 1 15 
» ar 4059) «0— Os 0-p- LES) Io 0-8) 19 


Dyz ta é) - (—10) 4 ta 2-04 70+) -0-F(148)-(-10) =5 


Il 


"ON, 1 1 1 1 7 
Sos ta (1-5) 04. 1-¢) 0+ 70 +e) 15— (L0) 1550 


Y ta 1) - (—10) Ta +n) -04 i000 50 n) -(—10) =0 


(6.33) 


Using the above formula to compute the derivative of the shape function N3 and 
substituting the above results, we obtain 


Ox 1 
" ON, » ON, » ON, n ON, 
Oy ( (x2. 3s) DA On 2 (X, Aca) (xs) 
, ON, NON, n ON, \ ƏN 
( Xa ðn v) 0E ( Nos OE 2 On 


a ON, VON n ON, \ ôN 
-( Dia tn) Ot + bor EE 2 E 


E j 5. (en) 0-748) 
(F 5-0-0) 0-1 +n) + oat Fé) 
1 
-(-4 an EF 
z 2 n n) E n n (6.34) 
g (1*9 99 + 9 


Moreover, it is seen that the determinant of the Jacobian matrix is 7 = D. 
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By substituting the derivatives of shape function N3 and the material parameters 
into equation (6.18), the expression for the first coefficient in K}, of element 1 is 


E [ON4 ƏN ƏN; ðN; 
Los tas a 0-9 d 20/2 |aa 


" D a3 ui = Ier +n)°(1— v)4 OK ey (he 20/2 jdédn 
7 h EDEN — 0.6) S) ¢ n)?-0.7 (a) gpa] © aan 


~ L art — 0.6) GS) c n) ^0,7 4- OK 4. ofa] Tačan 


= 1923076.92 f 
9, 


[0.02917 . (1+)? +0.01875 - (1+ 2d dédn 


(6.35) 


The integrand in the above formula is a quadratic polynomial of € and y. Thus, 
at least two integral points should be used in each direction, and there are four 
integral points in the element. Numerical integration is used to solve the problem: 


| E [ON ƏN; ON; ƏN; 
o, d | Ox Ox * Oy Oy 


(1 — 2v) Z dQ 


= 192307692. | [0.02917 - (1 +n)? +0.01875 - (1 + | dédy 
Qe 


1y? i’ 
= 1923076.92 - | 0.02917 - [1——2] +0.01875. (1— —] | -1-1 


. L (1+ A) +0.0875. (1 -5) TESI (6.36) 
+ L (: = AE (14 x) l ESI 


1^ px 
+ |0.02917- | 1+ — |] —-0.01875- (1+ —] |-1-1 
= 491480.77 


In this way, the stiffness coefficient can be obtained by numerical integration, 
further, other stiffness coefficients can be obtained similarly by the same computa- 
tion process. Through the above steps, the values of stiffness matrices of each ele- 
ment can be computed. 

For the global assembly of two-dimensional elements, the node numbers are 
marked on the left and upper sides of the matrix, and the expanded form of the 
stiffness matrix is given using the element location vector: 
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Node 1 2 3 4 5 6 
number 1 [K 2o K5 K3; K » Uy fi 
2 Ki, Ki, Ki; Kj, TD) f 
3 32 31 K3; tK) K3, +K3,) K5, K, ùs | _ J fs ; 
4 |K Kj (Kis T Kiz) (Kj, Ki) Ki, Ki ùf |f. 
5 K;, K;, K3; K5, us fs 
6 Kj Kj Kj Kj U6 fo 


where u and f are the global displacement and load matrices, respectively. The loads 
in the element can be obtained by integrating the product of the body forces and the 
shape functions: 


fi= f N,bdQ. (6.38) 
Q. 


Similarly, the node numbers are marked on the left of the load vector, and the 
expanded load is given using the element location vector: 


Node 1 hi ; 
number 2 f» fi 
3 Jf; ths 
Re ae : : 6.39 
4 Jf, itfi 639) 
5 [fs fy 
6 fs 4 


6.4 Treatments on Boundary Conditions 


Here, we focus on the displacement and traction boundary conditions for repre- 
sentative two-dimensional elastic problems. Other issues can be treated similarly at 
the corresponding nodes. 


(1) Displacement boundary conditions 

Using the above finite element form, it is very simple to impose displacement 
boundary conditions, as the parameters are now physical values, that is, they 
satisfy 


îi = à» = îy = às = ty = TW =O. (6.40) 
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Thus, the stiffness equation (6.37) can be expressed as: 


Node 1 2 3 4 5 6 
number 1 |Z 0 0 0 D 0 
2/0 I 0 0 tà» 0 
3 [0 0 (K&-K5) 00 O] Jal _ Jfsl, (6.41) 
4 |0 0 0 I. 0 0 p 0 
5 0 0 I1 0 ths 0 
6 0 0 0 I UG 0 


where I is the identity matrix. 
To simplify, rows and columns in the stiffness matrix of known displacement at the 
corresponding nodes can be deleted, thus reducing the size of the stiffness matrix: 


(Kł; + K2,) à; = f;. (6.42) 


'The matrices and vectors involved are expanded, leading to 


982961.54 0 tia hs 
| 0 pios] { v3 \ D Ea \ (6.43) 


(2) Traction boundary conditions 


To solve the above matrix equation, we should know the load at node 3. As there is 
no body force at present, the load matrix is 


(£)-() os 


Furthermore, the force acting on the boundary is considered. Imposing a traction 
condition at node 3 only requires the modification 


fy > fy tt=—1 x 10°. (6.45) 


The stiffness equation (6.42) to be solved becomes the following form: 


98296154 0 wi 4 
| vimos] l "s j 7 { —1000000 \ (6.46) 


which is solved directly, and the following results can be obtained 


{ N j H { ic \ (6.47) 


The displacement functions for each element can be obtained by combining the 
shape functions and node displacements: 
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4 
(£n) = M Nm ài 
a=1 


(6.48a) 
= NON 04 N04 N-0 
=0 
Element 1: 
4 
9 (£g) = 3 Nm 9 
ne (6.48b) 
= N,-0-- NS -04- N; - (—0.6592) + N, -0 
= —0.1648(1+ (1+1) 
4 
(En) — M ^ NalEn) à 
s (6.492) 
= N,-0--N5-04- N3-04- N4-0 
=0 
Element 2: 
4 
PE n) = Y Nn) 9 
i (6.49b) 


= N,-0+ Ny - (—0.6592) + N3 -0+ N4 -0 
= —0.1648(1 — &(1 +n) 


Once the displacements are known, the element stresses can be approximated as: 
4 4 
ô= M DS(N,) üj = X DB. (6.50) 
b=1 


In the computation process, the derivatives of each shape function are required. 
They can be computed as in the case of N3 above, as shown in equation (6.33). The 
derivatives of the other shape functions are 


ON, 1 ON» 1 
Oe E sd -g9 +0) 
am (=) 1 '" Jam (7) 1 | 
Oy 300-9) Vay} Vs? 
(6.51) 
aM) fia 
Ox o 30 | - n) 
ON, f — 1 
pm -3;048 
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The stress functions for each element can be obtained by combining the 
derivatives of the shape functions and node displacements: 
Element 1: 


ô! = [oi o, o1 Gene = = psi N,) à -Mona 


= D(Biüj + Bott, + B30; 2 um = D(B10 + B20 + B5ü; +B40) = DB; ix, 


ON; 
—- 0 
(1— v) U U 0 Ox 
_# s (1— v) v 0 0 A iu 
0 0  (1-2v)/2] | an, ƏN 
Oy Ox 
E 1- 0 os 0 
— RR 0 gj of ) 
v (=v) 0 0 j 0.6592 
0 0  (1=2v)/2]| 1. | 1 
55 Fé) 39 n) 
1 1 
39 ^ m - 9 ap (tt 9v 
1 1 
E sg * n? 99 0 * 90 — v) 0 | 
zt : i 
d tatn tarav 0.6592 
1 1 
gg (1. * 90 -20/2. 3c (1c (1 — 20)/2 
1 
gg * 9v 
1 
(1+ Q1 - v) 
= 965952 20 ! 
: L àv 
20 | 
tm — 20)/2 
0.15(1 + ¢) —19015.38(1+ €) 
EN OT! B ED SO) 
0.15(1 + ¢) —19015.38(1 + €) 
0.067(1 +n) —8451.28(1 +1) 


(6.52) 
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Element 2: 


P= È o& A J= Ypsi) a = pai 


= D(B à; + Bott; + B3 tt; EB. = D(B0 + Bott, + B30 +B,0) = DB», 


ON o 
(1— v) U U 0 Ox 
_E U (1— v) U 0 0 m as 
EE » (-9 0 5 ow 9 
0 0 0  (1-2v/21|80N, AN, 
Oy Ox 
Li 0 
(1— v) v U 0 30! +n) 
E E 0 Az 
== v — (1-v) v 0 350-8 { 
v (1— v) 0 0 0 
0 0 0 1 — 2v)/2 1 1 
S| qus a 
la 1 tii 
-ggQ + ma - ») 39 0-7 9v 
1 1 
E e UT) 39 4 QU - v) { 0 } 
"a 1 1 = 
d datas laor 0.6592 
1 1 
5504-90 -20/2 -zg 0 mü - 29/2 
(év 
i 0.15(1—£) 
z= (=é) — v) Z 
06592] 20 ” = 1267693314. 099019) 
d 1 0.15(1—£) 
a; (178v 
20 —0.067(1 + n) 
-390 * m = 29/2 
—19015.38(1—€) 
© J —44369.23(1-2) 
— | -19015.38(1—£) 
8451.28(1 +7) 
(6.53) 


Through the solutions of the above stress functions, the stress values at the 
integration points can be obtained by substituting the coordinates of the integration 
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points of the specified element. For example, the stress at the integral point (¢, n) = 


(- 44) of element 1 is 
19015 sa : ) 
l v3 
—19015.38(1 + é) MER ( 1 ) 
"mcm ; A 
a x y z ry = oa 4L S 
PAUSE UI) -19015.38(1 -) 
—8451.28(1 + 1) V3 
1 
—8451.28| 1 + — 
eras) 
—8036.84 
|] —18752.63 
—8036.84 
—13330.63 
(6.54) 


Furthermore, by substituting the node coordinates of the specified element into 
the solution of the above stress function, the stress values at the nodes can be 
obtained. For example, the stresses at node 2 of (é, n) = (—1, 1) of element 1 are as 
follows: 


—19015.38(1 + £) —19015.38(1 — 1) 
Mage Bt gh gba —44369.23(1+€) | — J —44369.23(1 — 1) 
Ty UR PN CES —19015.38(1 + £) —19015.38(1 — 1) 
—8451.28(1 +1) —8451.28(1 + 1) 
0 
7 0 
E 0 
—16902.56 
(6.55) 


6.5 Exercises 


6.5.1 Derive the weak form for general elasticity problems by following the steps 
below: 
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(a) Multiply each equation by an appropriate arbitrary function. 

(b) Integrate this product over the space domain of the problem. 

(c) Use integration by parts to reduce the order of the derivatives to a minimum. 
(d) Introduce boundary conditions if possible. 


6.5.2 In example 6.1, compute the value of the last coefficient in the stiffness matrix 
Ki, of element 1: 


E [ƏN; ON; ON; ðN; 
LUE Fe + ses 20/2| jan. (6.56) 


Chapter 7 
Solutions of Linear Algebraic Equations 


By using the finite element method, a standard discrete system is obtained. The 
basic stiffness equations of this system are linear. In general, a large number of such 
equations is difficult to solve by directly inverting the stiffness matrix. Herein, we 
introduce the commonly used LU decomposition method or triangular decomposi- 
tion method, also known as the Gaussian elimination method, including the forward 
elimination and back substitution procedures. 


7.1 LU Decomposition Method 


Consider a set of linear algebraic equations representing the global stiffness equa- 
tions in the finite element method, given by: 


Ki =f, (7.1) 


where K is an n X n square matrix with known values, representing the global 
stiffness matrix, & is a vector of unknown parameters, representing the global 
displacement vector, and f is a vector with known values, representing the global 
load vector. 

The matrix K is symmetric and positive definite; thus, it can be written as the 
product of a lower triangular matrix L with unit diagonals and an upper triangular 
matrix U as follows: 


K — LU, (7.2) 
where 
1 0 0 
M i i : (7.3) 
Lm Lao 1 
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and 
Ui, Ui» Uin 
0 Uz Von 
U= (7.4) 
0 0 PES Un 


This form is called a LU decomposition or triangular decomposition of K. The 
solution to the equations can now be obtained by sequentially solving the pair of 
equations (Demmel, 1997; Taylor, 1985; Strang, 1976; Wilkinson and Reinsch, 1971; 
Ralston, 1965): 


Ly=f (7.5) 
and 


Uü =y, (7.6) 


where y is introduced to facilitate the separation. In terms of the individual 
equations, the solution is given by: 


Yi = fis t=1 


il l (7.7) 
y =fi— Y Lb; i= 2, By ecg 
j=l 


and 


T 1= n 
ME C NM (7.8) 
sg (n- D w), i—-n-—1,n—2,..1 


where equation (7.7) is commonly called forward elimination, and equation (7.8) is 
called back substitution. 

As indicated in figure 7.1, it is convenient to subdivide the coefficient array in the 
stiffness matrix K into three parts: 


(1) Reduced zone: the region that is fully reduced. 

(2) Active zone: the region currently being reduced, where the jth column above 
the diagonal and the jth row to the left of the diagonal constitute the active 
zone. 

(3) Unreduced zone: the region that contains the original unreduced coefficients. 
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U  jth column active zone 


N 
Reduced 
zone 


jth row 
active zone 


Fic. 7.1 — Reduced, active, and unreduced zones in LU decomposition of stiffness matrix K. 


The coefficients in the decomposed matrices L and U can be stored in the active 
zone, respectively, as shown in figure 7.2. Through the above processes, the stiffness 
matrix is decomposed into two matrix parts. 


Fic. 7.2 — Matrices L and U stored in active zone. 


The algorithm for the LU decomposition of an n X n square matrix can be 
described using figure 7.3 as follows: 


104 
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Step 1. Active zone.First row and column to principal diagonal. 


Active zone 
Ky Ky Ky L=1 1 


' 
l 
psum deum 4 


i 
1 
Ky Ky Ky, : 
Kı Ka ^ Ks | 


Step 2. Active zone. Secondrow and column to principal diagonal.Use first row of K 
to eliminate L,,U,,: The active zone uses only values of K from the active zone and 
values of L and U which have already been computed in steps 1 and 2. 


y Reduced zone 
Active zone 
y 


Ky ! Kj à did 0 ] Ui | UK | 
Ka d Kat Ka [L-KEU. Ll [ _0 UeKeLUs: 
Ky Ky 1 


Step 3. Active zone. Thirdrow and column to principal diagonal. Use first row to 
eliminate L,,U,,: Use second row of reduced terms to eliminate L,,U,,(reduced 
coefficient K,,).Reduce column 3 to reflect eliminations below diagonal. 


Reduced zone 
a Active zone 
Y^ 


UyjzKy | 
Uy3=K3—L U3 i 
l 


i es a fas SS ss, aa aa a ee tee 


L,=K/U, 
Lyz(Ks- La Ui) Us 


Fic. 7.3 — LU decomposition of stiffness matrix K. 


The first row and column are activated as: 


Ui = Ku, Lu-l. 


For each active zone j from 2 to n, 


and 


K: 
La = Ds hy = Kij 
1 il 
Lii =F) Kj = Lim Umi 
1 m=1 
i-l 
Up = Ky- So beUa do339 uf 


and finally, with Lj; = 1, 


(7.11) 
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j-1 
Ujj = Kj = »» Lim Umi- (7.12) 


m=1 


Example 7.1 Using the triangular decomposition of the stiffness matrix K, solve the 
displacements in the stiffness equations: 


4 2 1 4 
Kü—f with K—|2 4 2 f=i2 (7.13) 
124 i 


First, the 3 X 3 stiffness matrix K is decomposed as shown in table 7.1. 
TAB. 7.1 — LU decomposition of 3 X 3 stiffness matrix K. 


K L U 
Stepl. Lii; Uj, = 4 
f 2 ] | | 4 | 
2.4 2 
124 
1 4 2 
L i | | | 


2 
Step2. Loi = i = 0.5; Uis = 2; Lo» =4 0.5x2—3 
2—0.25x2 21.5 


1 

2 2 

1 4 
3 3 


1 
Step3. La = i = 0.25; Uis = Ts La» = 
U23 —2—05x1- 1.5; L33 = 1; U33 =4-0.25x1-05x15=3 


1 1 
2 0.5 1 
12 4 0.25 0.5 1 


Step4. Check 
1 4 2 1 4 2 1 
0.5 1 3 15|=]|2 4 2 
0.25 0.5 1 3 1 2 4 


Second, using the decomposed matrices L and U, forward elimination and back 
substitution procedures are performed: 


t=1, y-fi-h-4 (7.14a) 


2 
4 
2 


= 0.5 


i-1l 1 
i=2, y=f- > Lyjy=h-YX Lan=2-05x4=0 (7.14b) 
jl jl 
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i—1 2 
i=3, yw=f—- do Ljy=h- X Lyy=1-0.25x4-05x0=0 — (7.140) 


j=l j=l 


and 
za Ea 4 
=3 ,=—=—=-=0 7.15 
"IU MEC Cg vn?) 
j= 2, (i 1 Y^ Us : Y` Unt (o 1.5 x 0) =0 
— ui = => pom ij Uy = — — 3 UZ = — — " — 
: ii Y ee: ve Uz2 9 = inn 3 
j=i+1 j=3 
(7.15b) 
"P Tons 1 ST 
i=, di — T com Uijit; "ES n - 5, Ui 
j=i+1 j=2 (7.15c) 
1 1 
E x Uiz ô + Uj3üg) -40-2 x0+1x0)=1 
The displacement solutions are obtained as follows: 
#={% à i} ={1 0 of. (7.16) 
As verified below, these solutions are completely reliable: 
4 2 1 1 4 
Ku=|2 4 2 O>=4 2>=f. (7.17) 
12 4 0 1 


Similar to the LU decomposition method introduced above, direct methods have 
been developed for solving the stiffness equations, such as LDLT and Cholesky 
decomposition methods. Users can use and implement any method according to 
their needs. 


7.2 Exercises 
7.2.1 Using the LU decomposition of the stiffness matrix K, compute the 


displacements: 


6 2 1 
Ki=f with K=|2 6 2 f=ias. (7.18) 
12 6 1 


Chapter 8 
Error Estimation and Adaptive Analysis 


The accuracy and effectiveness of finite element solutions depend on mesh quality. 
A few high-performance adaptive algorithms have been proposed to automatically 
optimise the mesh based on the complexity of the problem (Azadi and Khoei, 2011; 
Oden and Prudhomme, 2001; Zienkiewicz and Zhu, 1992a, 1992b; Babuška and 
Rheinboldt, 1978). Compared with the traditional method, the adaptive finite ele- 
ment algorithms exhibit improved accuracy and efficiency. It should be noted that 
adaptive algorithms have some obvious advantages over traditional numerical 
methods in terms of solving challenging problems, such as fracture, singularity, and 
eigenvalue problems. The adaptive finite element method uses the techniques of 
error estimation and refinement. 


8.1 Error Estimation of Finite Element Solutions 


8.1.1 Error of Finite Element Solutions 


'To improve the quality of finite element solutions, it is necessary to estimate the 
error involved. The finite element solution under the current mesh has the following 
error estimation: 


e — u—u, (8.1) 


where u^ and u are the conventional and exact solutions, respectively, and e is the 
error between them. However, it is difficult to obtain exact solutions for general 
problems. An approximate error estimate (a posteriori error estimator) can be 
obtained by replacing exact solutions with a superconvergent solution with higher 
accuracy than the conventional solution: 


e* = u*—ul, (8.2) 


where u* is the superconvergent solution, and e is the error between u* and the 
conventional solution. 
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8.1.2 Superconvergent Patch Recovery Method 


A typical error estimation method is based on superconvergent patch recovery 
(Zienkiewicz and Zhu, 1992a, 1992b), whereby it was demonstrated that the 
high-precision superconvergent points of stress solutions are Gauss integration points, 
and those of displacement solutions are element nodes. Figure 8.1 shows that the 
stress recovery of an element using high-precision stress results yields Gaussian inte- 
gration points as the natural superconvergent points in quadrilateral elements. The 
superconvergent stress solutions are obtained from the statistics of the patch nodes 
and superconvergent points of elements adjacent to the nodes, which can be used as an 


error estimator. 
2 |Patch recovery in 
adjacent elements 


Element: 1,2,3,4; Gaussian integration point: O 
Node: Patch node: @ 


Fic. 8.1 — Superconvergent patch recovery for nodal stresses of element. 


Furthermore, the superconvergent patch recovery displacement method was 
developed (Wang, 2021a; Wang et al., 2018; Wiberg et al., 1999a, 1999b) to acquire 
superconvergent finite element displacement solutions in static and dynamic prob- 
lems. As shown in figure 8.2, if element eis a superconvergent computation element 
and elements e — 1 and e + 1 can be considered its neighbouring elements, all finite 
element nodes in patched elements e — 1, e, and e + 1 are selected for the compu- 
tation process. The superconvergent displacements for element e can be computed as 


u'(z) = »IE Liam, (8.3) 


where r (=2) is the number of end nodes, s is the number of internal nodes and N; is 
the shape function matrix. Using the high-order shape function interpolation, the 
shape function polynomial order is increased, r+ s > m + 1. 


Fic. 8.2 — Superconvergent patch recovery for displacements of element e. 
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To optimise the superconvergent order O(h?™) for displacements at the end 
nodes, the displacement recovery field can be expressed for finite element nodes as: 


Ui(z)=Pa, i=1, ..., nq, (8.4) 


where P is the given function vector and a can be obtained by the least-squares 
fitting technique for the coincidence of displacements at the end nodes in both the 
recovery and the conventional finite element fields. The superconvergent displace- 
ments of the recovery field are used in equation (8.3) to obtain the superconvergent 
solutions on element e. The vector coefficients P and a forms used for this study 
were: 


P=|1 £u æ], a-—[a wd oe am], (8.5) 


The value of a was determined by obtaining the minimum value of the following 
functional so that the product result on the positions of the finite element nodes in 
equation (8.5) is equal to the displacement values: 

n 
II = (u(zs)-P(z)a) i=1, ..., na, (8.6) 


FA 


where n is the node number of all elements patched together. 
Using the least square method to solve equation (8.7), the value of the coefficient 
a can be obtained as: 


a — A lb, (8.7) 


The coefficient matrices A and b are defined as follows: 
A=) P(u)P, b=> P(x) um. (21,2 na. (8.8) 
j=l j=l 


After the coefficient a is determined, the superconvergent solutions of the dis- 
placement of the piecewise elements can be obtained from equation (8.3). 

To determine if the solution on the current mesh satisfy the given tolerance, the 
error is estimated by 


: 2 g 1/2 
lel < Tot- [CI]. + hel?) /n] (8.9) 
where n, is the number of elements, and |le*|| is an error in energy norm: 
1/2 
lle*l| = [a(€*, e] ?— |/ eT Le i (8.10) 
Q 


where e* = u*—u", and Q is the solution domain. Equation (8.9) can be equivalently 
rewritten in the following form: 


é= with e= Tol: (fal? het?) /.] . (8.11) 
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where € should satisfy: 
ex]. (8.12) 


8.2 Adaptive Finite Element Method 


8.2.1 Categories of Adaptive Finite Element Method 


The adaptive finite element method can be implemented by the following three-step 
adaptive procedures: 


(1) FE solutions. On the current mesh, the finite element solution u^ is obtained by 
the conventional finite element method. 

(2) Error estimation. Based on the finite element solution, the superconvergent 
solution u* is computed through superconvergent patch recovery displacement. 
Subsequently, u* is used to estimate the error of the u^, and then the errors in 
energy norm on the global domain are marked out. 

(3) Solution refinement. If the pre-specified error tolerance is not satisfied for the 
elements above, the mesh is further subdivided to generate a new mesh based on 
the computed errors, and the mesh refinement procedure is used to refine the 
current mesh on the global domain. Then, the procedure returns to the first step 
(i.e., FE solutions), and the cycle repeats until all the elements satisfy the 
pre-specified error tolerance. 


Solution satisfying 
pre-set Tolerance 


Satisfies 
Tolerance 


Solution refinement 
Does not 


satisfy Tolerance 


Fic. 8.3 — Basic process flow of the high-performance adaptive finite element algorithms. 


There are various procedures for the refinement of finite element solutions, and 
they broadly fall into three categories: 


(1) h-version refinement: The same order of element polynomial is used, but the 
elements are changed in size (in some locations, they are made larger, and in 
others, smaller) to achieve maximum efficiency in reaching the desired solution. 

(2) p-version refinement: The procedure uses the same element size and simply 
increases the order of the element polynomials. 

(3) hp-version refinement: The procedure combines the h- and p-version refine- 
ments. In this procedure, the element size and the degree of the element poly- 
nomial are altered. 
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The implementation of h-refinement is relatively simple. Based on the error 
estimation of the superconvergent solution introduced above, a developed h-version 
adaptive finite element method is given below. 


8.2.2 h-Version Adaptive Finite Element Method 


If equation (8.12) is not satisfied, the corresponding element needs to be subdivided 
into uniform sub-elements by inserting some interior nodes through h-refinement. 
These are computed by: 


Inew = EVP) ia; (8.13) 


where Anew is the length of the sub-element and A is the original length of element 
e. The above element subdivision approach is implemented as follows: 


= min( Em d), (8.14) 


where Mow is the number of sub-elements after element subdivision, the symbol |-| 
represents the ‘floor’ operator (i.e., the rounding up to the nearest integer), and d is 
the limit number needed to avoid too many redundant elements. Each element 
e that does not satisfy the pre-specified error tolerance is uniformly subdivided by 
h-version mesh refinement, e.g., hnew = haa/6, as shown in figure 8.4. 


Pow 
j ju @ Initial FE node 
e 
le | ® Now FE node 
has 


Fic. 8.4 — Diagram of uniform subdivision on element e by h-version mesh refinement (e.g. 


liey = haa/6). 


The use of non-uniform meshes in the high-performance adaptive analysis is a 
crucial factor for efficient computation, where a minimum number of elements and 
optimised node distribution of nodes are used to quickly reduce errors and obtain 
high-precision solutions. Based on the error estimation in the energy norm, the solu- 
tion error of an element can be estimated using the superconvergent solution for each 
element, and it can be effectively reduced by using a uniform subdivision for this 
element. Moreover, to avoid redundant elements in this uniform subdivision, the 
adaptive finite element method limits the maximum number of subdivisions in each 
adaptive stage. Overall, the mesh in the global domain is non-uniform because each 
element is independent of the error estimation and mesh subdivision. With this sub- 
division and the element increase, the non-uniform tendency of the elements becomes 
prominent to suit displacement variation. Through numerical examples in which 
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disturbed displacements around cracks and highly varying displacements are obtained 
by adaptive analysis with optimised meshes, as shown in Figures 8.5 and 8.6 
(Wang, 2021b; Wang et al., 2018), it is demonstrated that the error estimation in the 
energy norm and the uniform subdivision for each element for mesh refinement allows 
for high accuracy, reliability, and effectiveness in obtaining displacement solutions. 


Fic. 8.5 — Computed solutions of disturbed displacements around cracks by adaptive analysis, 
and the final adaptive meshes are shown as tick marks on the horizontal axis. 


10.0 
5.0 

= 0.0 ul, 
-5.0 
-10.0 


[pr pr npn np JANJAN LE G4 — 
0.0 02 04 0.6 0.8 1.0 
x/l 


Fic. 8.6 — Computed solutions of dramatically changed displacements by adaptive analysis, 
and the final adaptive meshes are shown as tick marks on the horizontal axis. 


8.2.3 hp-Version Adaptive Finite Element Method 


Based on the h-version adaptive FEM introduced above, a novel Ap-version adaptive 
FEM can be developed (Wang and Wang, 2022). This section introduces the basic 
algorithm and computation procedure of this hp-version adaptive FEM and shows 
the effectiveness of the Ap-version refinement and the reliability of the results. The 
hp-version adaptive finite element algorithm is proposed for determining the 
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solutions of shell structures via error homogenisation and higher-order interpolation. 
This algorithm first develops the established h-version mesh refinement method for 
detecting the non-uniform distributed optimised meshes, where the error estimation 
and element subdivision approaches based on the superconvergent patch recovery 
displacement method are introduced to obtain high-precision solutions. The errors 
in the displacement solutions in the global space domain are homogenised and 
approximately the same. Subsequently, on the refined meshes, the algorithm uses 
higher-order shape functions for the interpolation of trial displacement functions to 
reduce the errors quickly until the solution meets a pre-specified error tolerance 
condition. In this algorithm, the non-uniform mesh generation and higher-order 
interpolation of shape functions are suitable for addressing the problem of complex 
displacements caused by variable structural geometries. 

This Ap-version refinement procedure involves the following four basic 
techniques: 


(1) FE solution: Conventional FE computation is performed based on the initial 
mesh. The FE solution u” is obtained under the current mesh. 

(2) Error estimation: Based on the FE solution w^, the superconvergent solution u* 
is obtained using the superconvergent algorithm of FE post-processing. At the 
same time, the error estimates between the superconvergent and FE solutions 
can be obtained. 

(3) h-version mesh refinement using error estimation and element subdivision: For 
the element whose error estimation does not meet the error tolerance condition, 
the mesh subdivided refinement method is used to subdivide it to obtain a new 
FE mesh before returning to step (1). If all the elements meet the error tolerance 
condition, no subdivision is necessary, and the solution process is completed. 

(4) p-version enrichment using higher-order shape functions for interpolation: 
Under the optimised mesh obtained in the above steps, the error tolerance is 
reduced, the element order is increased, the FE computation continues, and the 
new error tolerance is used to estimate the error. For the element whose error 
estimate does not meet the requirements, this step is repeated, and the order of 
the element is increased until the new specified error tolerance condition is met, 
completing the solution process. 


To verify the efficiency of the hp-version adaptive FEM, the error distributions of 
the computed solutions of displacements for a representative shell structure in mesh 
refinement and higher-order interpolation stages have been analysed. Figure 8.7 
shows the error distributions of the computed displacements in the mesh refinement 
and precision control stages of the hp-version adaptive FEM. To facilitate com- 
parison, the error results of corresponding mode components are listed in 
Figure 8.7a and b, respectively. For a pre-specified error tolerance level Tol = 1079, 
the computed displacement in the h-version method exceeds the requirement. The 
refined meshes in the whole domain yield a solution error distribution of displace- 
ment that is almost a running average of the tolerance Tol,. After the order of the 
shape function of the element is increased by 1, all the modes of the hp-version 
method strictly meet the error limitation requirements. Thus, the Ap-version 
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refinement may use fewer optimised meshes than the h-version mesh refinement, and 
only one-step interpolation of the higher-order shape function can provide the dis- 
placement satisfying the accuracy requirement. 


-7 
x10° 8.0)x 10 
8.0 
4.0 4.0 
* 0.0 w 0.0 
-4.0 40 
-8.0 : 
0.0 0.2 0.4 0.6 0.8 1.0 gaan ILA J- 
vi 0.0 0.2 0.4 0.6 0.8 1.0 
x/l 
(a) Error of computed solutions of displacements (b) Error of computed solutions of displacements after 
after mesh refinement, and the adaptive meshes higher-order interpolation, and the adaptive meshes 
are shown as tick marks on the horizontal axis. are shown as tick marks on the horizontal axis. 


Fic. 8.7 — Error distributions of the computed solutions of displacements in mesh refinement 
and higher-order interpolation stages. 


8.3 Exercises 


8.3.1 Summarise the basic methods and procedures in adaptive analysis. 


Chapter 9 
Programs of Finite Element Method 


Based on the basic methods and theories introduced in previous chapters, in this 
chapter, we introduce finite element computer programs and discuss the input 
parameters and the output for each example. It should be noted that we do not focus 
on the specific programming implementation, but along with the source code, we 
also directly provide the executables, thus aiding readers in understanding and using 
the techniques. Readers who have an in-depth understanding of specific program- 
ming methods and more source codes may refer to the related literature devoted to 
the programming and programs of the finite element method (Smith et al., 2013; 
Zienkiewicz et al., 2013b). The programs described in this chapter are written in 
FORTRAN, and the compilers used include Intel Visual Fortran, G95, GFOR- 
TRAN, Compaq Visual Fortran, Cray, NAG, and the Portland Group. Readers may 
use the compiler of their choice. 

We provide examples of one-, two-, and three-dimensional finite element pro- 
grams, including the one-dimensional program for beam deformation, the 
two-dimensional program for plane strain problems, and the three-dimensional 
program for solid compression. The programs are provided in the following forms: 
(1) As source code so that readers can understand the logical relationship of the 
program, further improve the code, and realise new functions; for illustration, the 
main section of each program is presented and described in detail below. (2) As 
compiled executable files so that readers can ignore the specific implementation, run 
the program directly, and obtain the results. 


9.1 One-Dimensional Program of Beam Deformation 


This program is suitable for example 3.2, and it uses the geometric model of the 
one-dimensional program of beam deformation shown in figure 9.1. It solves the 
cantilever beam problem in that example, with Young's modulus E = 1000, length 
a = 10, loading t, = 25, and loading b, = 10 (0 < x < a/2). 
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b=10 
fos o Ó o ——————- 
1 © 2 Q 3 Q9 4 ©) 5 
a a a a 
| »| 
Element number: @~@ Node number: 1~5 


Fic. 9.1 — Geometric model of the one-dimensional program of beam deformation. 


9.1.1 Main Program 


The main program of the one-dimensional program of beam deformation is shown 
below and contains the following modules: dynamic arrays, input, and initialisation, 
loop the elements to find global arrays sizes, global stiffness matrix assembly, read 
loads and/or displacements, equation solution, and retrieve element end actions. 


PROGRAM p1DFEM 

USE main 

USE geom 

IMPLICIT NONE 

INTEGER, PARAMETER:: ivp-SELECTED REAL KIND(15) 

INTEGER:: fixed freedoms, i, iel, k, loaded nodes, ndof=2, nels, neq, nlen, nod=2, & 
nodof=1, nn, nprops=1, np types, nr 

REAL(iwp):: penalty-1.0e20 iwp, zero—-0.0 iwp 

CHARACTER(LEN-15):: argv 

]l---------- dynamic arrays -————-—-—--------- 

INTEGER, ALLOCATABLE:: etype(:), g(:), g g(: :), kdiag(:), nf(:, :), no(:), & 
node(:), num(:) 

REAL(iwp), ALLOCATABLE:: action(:), eld(:), ell(:), km(:, :), kv(:), loads(:), & 
prop(:, :), value(:) 

]l--------- input and initialisation -——————-—-—--—-— 

CALL getname(argv, nlen) 

OPEN(10, FILE=argv(1:nlen) //’.dat’) 

OPEN(11, FILE=argv(1:nlen)//’.res’) 

READ(10, *) nels, np_ types 

nn=nels+1 

ALLOCATE(g(ndof), num(nod), nf(nodof, nn), etype(nels), ell(nels), eld(ndof), & 
km(ndof, ndof), action(ndof), g g(ndof, nels), prop(nprops, np types)) 

READ(10, *) prop 

etype=1 

IF(np_types>1) READ(10, *) etype 

READ(10, *) el 

nf=1 

READ(10, *) nr, (k, nf(:, k), i=1, nr) 

CALL formnf(nf) 

neq=MAXVAL (nf) 
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ALLOCATE(kdiag(neq), loads(0:neq)) 
kdiag=0 
|o—e-e—e———- loop the elements to find global arrays sizes — — 
elements 1: DO iel=1, nels 
num=(/iel, iel+1/) 
CALL num_to_g(num, nf, g) 
g_g(:, iel)=g 
CALL fkdiag(kdiag, g) 
END DO elements_ 1 
DO i=2, neq 
kdiag(i)=kdiag(i)+kdiag(i-1) 
END DO 
ALLOCATE(kv(kdiag(neq))) 
WRITE(11, ‘(2(A, I5))) & 
" There are", neq, " equations and the skyline storage is", kdiag(neq) 
| global stiffness matrix assembly- - — — — — — — 
kv=zero 
elements _2: DO iel=1, nels 
CALL rod km(km, prop(1, etype(iel)), ell(iel)) 
&—g 8g; iel) 
CALL fsparv(kv, km, g, kdiag) 
END DO elements 2 
! read loads and/or displacements - - ——————— 
loads=zero 
READ(10, *) loaded_nodes, (k, loads(nf(:, k)), i=1, loaded nodes) 
READ(10, *) fixed freedoms 
IF (fixed freedoms/—0)THEN 
ALLOCATE(node(fixed freedoms), no(fixed_ freedoms), value(fixed freedoms)) 
READ(10, *) (node(i), value(i), i=1, fixed_ freedoms) 
DO i=1, fixed_ freedoms 
no(i)=nf(1, node(i)) 
END DO 
kv(kdiag(no))=kv(kdiag(no))+-penalty 
loads(no)=kv(kdiag(no))*value 
END IF 
l= Sasa ase esas cquationmisolution m mme 
CALL sparin(kv, kdiag) 
CALL spabac(kv, loads, kdiag) 
loads(0)=zero 
WRITE(11, ‘(/A)’)" Node Disp" 
DO k=1, nn 
WRITE(11, (I5, 2E12.4)’) k, loads(nf(:, k)) 
END DO 


lm reirievetelemenutend(actlorse E 
WRITE(11, ‘(/A)’) " Element Actions" 
elements 3: DO iel=1, nels 
CALL rod km(km, prop(1, etype(iel)), ell(iel)) 
&—g g(; iel) 
eld=loads(g) 
action2MATMUL(km, eld) 
WRITE(11, ‘(15, 2E12.4)’) iel, action 
END DO elements 3 
STOP 
END PROGRAM p1DFEM 
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9.1.2 Numerical Example 


The input parameters for the above geometric model and the basic conditions are 
contained in file IDFEM.dat; four elements with uniform length are used; the input 
data and their interpretation are listed in table 9.1. 


TAB. 9.1 — Input data and interpretation of the one-dimensional program of beam 
deformation. 

Input data Interpretation 

4 1 element in x, material 
1000.0 EA 

2.5 25 2.5 2.5 element length 

1 nr displacement BC 
10 k, nf(:,k), i—1, nr 

5 nr traction BC 

1125 225 3125 40 525 k, nf(:,k), i=1, nr 

0 fixed freedoms 


Based on running the one-dimensional program of beam deformation, the output 
solutions of the program for the above geometric model and basic conditions are 
contained in file IDFEM.res, and are shown below: 


Node Disp 
1 0.0000E-4-00 
2 0.1563E+00 
3 0.2500E4-00 
4 0.3125E+00 
5 0.3750E4-00 

Element Actions 
1 -0.6250E+02 0.6250E+02 
2 -0.3750E+02 0.3750E--02 
3 -0.2500E+02 0.2500E--02 
4 -0.2500E4-02 0.2500E+02 


9.1.3 Interactive Interface 


To facilitate learners using the one-dimensional (1D) program of beam deformation 
for analysis, this book provides a simple interactive interface for learners to input 
parameters, compute, and view results. This interface is displayed and explained in 
sequence according to the following steps: 

First, the interface is opened (figure 9.2), which exhibits the message ‘Select the 
program to run’. The option ‘One-dimensional program of beam deformation’ is 
selected for analysing the 1D problem. Similarly, two-dimensional (2D) and 
three-dimensional (3D) problem types can also be selected by selecting the corre- 
sponding options. 
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fy. AFEMRS 1.5 


m 


Select the program to run: 


One-dimensional program of beam deformation | 


Two-dimensional program of plane strain problem 


Three-dimensional program of solid compression 


Dr. Yongliang Wang, www.wangyongliang net, wangyl@cumtb.edu.cn 


Fic. 9.2 — Interface of selection for the one-, two-, and three-dimensional programs. 


Figure 9.3 shows that this interface includes the parameter input window on the 
upper left, the model diagram on the upper right, the summarised input data on the 
lower left, the output data in the lower middle part, and the displacement diagram 
on the lower right. 


E One-dimensional program of beam deformation - oF x 
Element number Material 
Tensile stiffness EA b 
Element length ^a | 
Dis. BC no. 
Node number of Value of Dis. BC aa 
Tra. BC no. EEEND a + b c - d 
Node number of Tra. BC. Value of Tra. BC Add 
Fei geste Element number: ©- Node number: 1-5 
Model of one-dimensional program of beam deformation 
Display input data | Generate fle IDFEM dat | Compute 
Display output data Plot Retum Exit 
[Display input data: Display output data: Displacement diagram: 


1.00 
0.90 
0.80 
0.70 
0.60 
0.50 
0.40 
0.30 
0.20 
0.10 


01234567 8 9 10 


Fic. 9.3 — Interface of the one-dimensional program of beam deformation. 
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(1) Step 1 


Figure 9.4 shows two parameters in the first line of the parameter input window: 

Element number: number of elements used in the z-direction. The value entered 
in this figure is 4. 

Material: number of material types. The value entered in this figure is 1. 


© One-dimensional program of beam deformation = 0 X 
Element number 4 Material 1 
Tensile stiffness EA b 
Add 
oe a 
D 102 0 3 9 4 @ 5 
Node number of Value of Dis. BC Add | | ] | 
Tra. BC no. s ii b N ii I 4 
Node number of Tra. BC Value of Tra. BC Eni) 
Fixed freedom Element number: ©~@ Node number: 1-5 
Model of one-dimensional program of beam deformation 
Display input data Generate file IDFEM.dat Compute 
Display output data Plot Retum Exit 
Display input data: Display output data: Displacement diagram: 
1.00 
0.90 
0.80 
0.70 
0.60 
0.50 
0.40 
0.30 
0.20 
0.10 
0123245267 8 9 10 
Fic. 9.4 — Element number and material in the one-dimensional program of beam 


deformation. 


(2) Step 2 


Figure 9.5 shows one parameter in the second line of the parameter input 
window: 

Tensile stiffness EA: tensional stiffness of the beam. The value entered in this 
figure was 1000. 


(3) Step 3 


Figure 9.6 shows one parameter in the third line of the parameter input window: 

Element length: length of each element. The length value of each element can be 
input sequentially, and the option *Add' is selected. Subsequently, the window *Add 
successfully’ appears, and the option ‘Agree’ is selected to confirm the addition. The 
values entered in this figure were 2.5, 2.5, 2.5, and 2.5. Notably, the model diagram 
shows that the length of each element is the same (a — b — c — d); regarding the 
general input parameters and computation, the length for different elements can be 
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IS One-dimensional program of beam deformation - n 


n —] 
Add 


Tensile stiffness £4 1000 


© 2 © 3 © © s 
Value of Dis. BC Add _ _| 
Tra. BC no. - A e q 


Node number of Tra. BC Value of Tra. BC Ada. 


TEE] 


Element number: (1)-(8) Node number: 1-5 


Fixed freedom 
Model of one-dimensional program of beam deformation 
Display input data | Generate le IDFEM dat | Compute. | 
Display input data: Display output data: Displacement diagram: 


0123 *5 67 8 9 10 


Fic. 9.5 — Tensile stiffness in the one-dimensional program of beam deformation. 


Tensile stiffness EA 1000 


Element length 25 


Dis. BC no. 
Node number of Value of Dis. BC Add 
Tra. BC no. 
Node number of Tra. BC Value of Tra. BC. ada 
Fixed freedom © One-dime.. — O — X jtnumber: (0-5 Node number: 1~5 
lof one-dimensional program of beam deformation 
Display input data | Generate file IDFEM dat | uS 
Display output data. Plot Return 
IDisplay input data: Display Displacement diagram: 
1.00 
0.90 
0.80 
0.70 
0.60 
0.50 
0.40 
0.30 
0.20 
0.10 
01235354567 8 9 10 


Fic. 9.6 — Element length in the one-dimensional program of beam deformation. 
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set to different values. Notably, after adding each group of parameters successfully, 
this line of parameters is reserved to facilitate the generation of more elements with 
the same length, which can be manually cleared to input the next group of 
parameters. 


(4) Step 4 


Figure 9.7 shows one parameter in the fourth line of the parameter input 
window: 

Dis. BC no.: number of displacement boundary conditions. The value entered in 
this figure is 1. 


I One-dimensional program of beam deformation - D x 


Element number 4 Material 1 
Tensile stiffness E4 1000 
Element length 2.5 Add 
Dis. BC no. : 
Node number of Value of Dis. BC ERR 
Tra. BC no. 
Node number of Tra. BC Value of Tra. BC 2n 
Fixed freedom Element number: ©~@ Node number: 1-5 
Model of one-dimensional program of beam deformation 
Display input data Generate file IDFEM. dat. Compute 
Display output data Plot Retum Exit 


Display input data: 


Display output data: 


Displacement diagram: 


1.00 
0.90 
0.80 
0.70 
0.60 
0.50 
0.40 
0.30 
0.20 
0.10 


01 z 3 43 67 8H 9 10 


Fic. 9.7 — Number of displacement boundary conditions in the one-dimensional program of 
beam deformation. 


(5) Step 5 


Figure 9.8 shows two parameters in the fifth line of the parameter input window: 

Node number of Dis. BC: node number of the displacement boundary, specifying 
the node at which a specific displacement occurred. The value entered in this figure 
is 1. 

Value of Dis. BC: the value of the displacement boundary, specifying that used 
for the displacement boundary. The value of the displacement boundary of each 
specifying node can be input sequentially, and the option *Add' is selected. Then, 
the window ‘Add successfully’ appears, and the option ‘Agree’ is selected to confirm 
the addition. The value entered in this figure was 0. 
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Node number: 1~S 


f one-dimensional program of beam deformation 


5 

Element number 4 Material 1 

Tensile stiffness EA 1000 

Element length 2.5 Add 

Dis. BC no. 1 

Node number of 1 Value of Dis. BC |0 | 

Tra. BC no. a —+\-— b 

Node nuniber of Tra. BC Value of Tra. BC Add 

Fixed freedom B One-dime.. — X jtnumber: ®©- 

Display input data Generate file IDFEM dat peru UE 
Display output data Plot Retum ——J 

Display input data: Display 
1.00 
0.90 
0.80 
0.70 
0.60 
0.50 
0.40 
0.30 
0.20 
0.10 


Displacement diagram: 


0123456 7 8 9 10 


Fic. 9.8 — Node number and value of the displacement boundary condition in the 
one-dimensional program of beam deformation. 


(6) Step 6 


Figure 9.9 shows one parameter in the sixth line of the parameter input window: 
Tra. BC no.: number of traction boundary conditions. The value entered in this 


figure is 5. 


1.00 
0.90 
0.80 
0.70 
0.60 
0.50 
0.40 
0.30 
0.20 
0.10 


E One-dimensional program of beam deformation = O X 
Element number 4 Material 1 
Tensile stiffness £4 1000 
Element length 25 Add 
Dis. BC no. 1 
Node number of 1 Value of Dis. BC | Add 
Tra. BC no. 5 
Node number of Tra. BC Value of Tra. BC Aaa 
Fixed freedom Element number: ©~@ Node number: 1-5 
Model of one-dimensional program of beam deformation 
Display input data Generate file IDFEM.dat Compute 
Display output data | Pot | Reus | Exi 
Display input data: Display output data: Displacement diagram: 


0123245267 8 9 10 


Fic. 9.9 — Number of traction boundary conditions in the one-dimensional program of beam 


deformation. 
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(7) Step 7 


Figure 9.10 shows two parameters in the seventh line of the parameter input 
window: 

Node number of Tra. BC: node number of the traction boundary, specifying the 
node at which specific traction occurs. The values entered in this figure were 1, 2, 3, 
4, and 5. 

Value of Tra. BC: value used for the traction boundary. The value of the traction 
boundary of each specifying node can be input sequentially, and the option ‘Add’ is 
selected. Subsequently, the window ‘Add successfully’ appears, and the option 
‘Agree’ is selected to confirm the addition. The values entered in this figure were 
12.5, 25, 12.5, 0, and 25. 


5 
Element number 4 Material 1 
Tensile stiffness E4 200 b 
Element length 25 Add 
Dis. BC no. i y 
Node number of 1 Value of Dis. BC | Add | 
Tra. BC no. 5 EM peesee sessi 
Node number of Tra. BC |1 Value of Tra BC|125 ^ | Add 
Fired keton D One-dime. — O X i number: -O Node number: 1-5 
= = of one-dimensional program of beam deformation 
Display input data Generate file IDFEM dat Add successfully 
Display output data Plot Retum 
Display input data: Display | Displacement diagram: 
] 1.00 
0.90 
0.80 
0.70 
0.60 
0.50 
0.40 
0.30 
0.20 
0.10 
D-t$t2g4ssvnuT 89 10 
Fic. 9.10 — Node number and value of the traction boundary condition in the 


one-dimensional program of beam deformation. 


(8) Step 8 


Figure 9.11 shows one parameter in the eighth line of the parameter input 
window: 

Fixed freedom: number of nodes with specified displacement. The value entered 
in this figure is 0 by default. 
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E3 One-dimensional program of beam deformation = p 
Element number 4 Material 1 
Tensile stiffness E4 1000 b, 
2.5 Add 
ae et , 
: 1 j, o 4 UNS 
MEC 102 © 5 © 4 O's 
Node number of i Value of Dis. BC | Ada 
Tra. BC no. 5 E : a 
Node number of Tra. BC |S Value of Tra. BC |25 eee 
Fixed freed 0 Element number: ©~@ Node number: 1-5 
Model of one-dimensional program of beam deformation 
Display input data | Generate fle IDFEM dat | Compute 
Display output data Plot Retum Exit 
Display input data: Display output data: Displacement diagram: 


1.00 
0.90 
0.80 
0.70 
0.60 
0.50 
0.40 
0.30 
0.20 
0.10 


012 345 6&7 8 9 10 


Fic. 9.11 — Fixed freedom in the one-dimensional program of beam deformation. 


(9) Step 9 


As shown in figure 9.12, the option ‘Display input data’ is selected, and the 
aforementioned input dataset is summarised on the lower left of the interface. 


5 
112.5 225 312.5 40 525 
lo 


& One-dimensional program of beam deformation - o x 
Element number 4 Material 1 
Tensile stiffness E4 anue b 
2.5 Add 
ose =a 
8 1 A o — 
Dn: 10252 0 3 9 4 @ 5 
Node number of 1 Value of Dis. BC |0 Add 
Tra. BC no. 5 a i 5 i Ñ g 
Node number of Tra. BC |5 Value of Tra. BC |25 aii 
Fixed freedom lo Element number: 0- Node number: 1-5 
Model of one-dimensional program of beam deformation 
Generate le IDFEM dat Compute 
Plot Retum Exit 
hai Display output data: Displacement diagram: 
1000 100 
5 25 25 25 e 
i 5 2.5 25 2. 090 
0.80 
10 0.70 


0.60 
0.50 
0.40 
0.30 
0.20 
0.10 


0123456789 10 


Fic. 9.12 — Display of input data in the one-dimensional program of beam deformation. 
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(10) Step 10 


As shown in figure 9.13, the option ‘Generate file IDFEM.dat' is selected, and 
the input dataset is recorded in the file 1DFEM.dat, which is saved in the path of the 
current file. Subsequently, the window ‘Generate successfully’ appears, and the 


option ‘Agree’ is selected for confirmation. 


1000 

25 25 25 25 
1 
10 


1.00 
0.90 
0.80 
0.70 
0.60 
0.50 


5 
Element number 4 Material 1 
Tensile stiffness EA 1000 
Element length 25 Add 
Dis. BC no. : 
Node number of 1 Value of Dis. BC |0 Add 
Tra. BC no. 5 
Node number of Tra. BC |5 Value of Tra. BC |25 Add | 
Fixed freedom o f One-dime.. — O — X jtnumber: ©- Node number: 1~5 
f one-dimensional program of beam deformation 
Display input data | Generate fle IDFEM dat Ea 
Display output data Piot | Renum | 
41 Display Displacement diagram: 


5 
1125 225 312.5 40 525 
jo 


0.40 
0.30 
0.20 
0.10 


6 3 2 3 4567 89 19 


Fic. 9.13 — Input data generation in the one-dimensional program of beam deformation. 


(11) Step 11 


As shown in figure 9.14, the option ‘Compute’ is selected, and the computation 
program starts to conduct operations using the aforementioned input data. Once 
the computation is completed, the window ‘Compute successfully’ appears, and the 
option ‘Agree’ is selected for confirmation. The output dataset above is recorded in 
the file 1DFEM.res, which was saved in the path of the current file. 


(12) Step 12 


As shown in figure 9.15, the option ‘Display output data’ is selected, and the 
output data of computed results are displayed on the lower right of the interface. 


Programs of Finite Element Method 127 


Bs EE EPE A prog of beers T TEREF- - 
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Tensile stiffness EA 1000 b 

Element length 2.5 Add T i í 
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Node number of Tra. BC |5 Value of Tra. BC |25 . Ada | 


Fixed freedom o D One-dime.. — O — X jtnumbe: O~@ Node number: 1~5 
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41 Display Displacement diagram: 
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Fic. 9.14 — Computation control in the one-dimensional program of beam deformation. 


I One-dimensional program of beam deformation - D x 
Element number s Material u 
Tensile stiffness EA 1000 k 


Dis. BC no. 1 


a mum 
jo o ran o - 
Node number of 1 Value of Dis. BC |0 Add | | | | 
a — b c + d 


Tra. BC no. d 
Node number of Tra. BC |5 Value of Tra. BC |25 e 
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Fic. 9.15 — Output data display in the one-dimensional program of beam deformation. 


(13) Step 13 


As shown in figure 9.16, the option ‘Plot’ is selected, and the diagram based on 
the displacement results is shown. Once the diagram is completed, the window 
‘Draw successfully’ appears, and the option ‘Agree’ is selected for confirmation. 
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A control option called ‘Return’ is also present in the interactive interface, which 
returns to the interface for selecting the computation program shown in figure 9.2. 
The control option ‘Exit’ exits the current interactive interface. 

The aforementioned operation steps are introduced only through a specific 
example. Users can use and analyse a variety of examples. Through the resulting 
data computed by this interface, users can also use other software to visually view 
the model and results. 


5 

Element number 4 Material 1 

Tensile stiffness EA 1000 

Element length 25 Add 

Dis. BC no. 1 
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Tra. BC no. 5 
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4 -0.2500E+02 0.2500E+02 


Fic. 9.16 — Result diagram in the one-dimensional program of beam deformation. 


9.2 Two-Dimensional Program of Plane Strain Problem 


This program is suitable for example 6.1. It uses the geometric model of the 
two-dimensional program of the plane strain problem shown in figure 9.17, with 
dimensions a = 30 in the horizontal x-direction and b = 10 in the vertical ydirec- 
tion. The material parameters are Young's modulus E = 1 X 10° and Poisson's ratio 
v — 0.3. The boundary conditions are that the upper left and upper right vertices 
are fixed, and the lower parts are fixed. The vertical concentrated load Facting on 
the middle domain is 1 x 10°. 


9.2.1 Main Program 


The main program of the two-dimensional program of plane strain problem is shown 
below and contains the following modules: dynamic arrays, input, and initialisation, 
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Element number: @ and @) Node number: 1~6 


Fic. 9.17 — Geometric model of the two-dimensional program of plane strain problem. 


loop the elements to find global arrays sizes, element stiffness integration and 
assembly, equation solution, recovery stresses at nip integrating points, and recover 
stresses at node points. 


! Program 2 Two-dimensional program of plane strain problem 

PROGRAM p2DFEM 

USE main 

USE geom 

IMPLICIT NONE 

INTEGER, PARAMETER:: ivp-SELECTED REAL KIND(15) 

INTEGER:: fixed freedoms, i, iel, k, loaded nodes, ndim=2, ndof, nels, neq, nip, & 
nlen, nn, nod, nodof=2, nprops=2, np types, nr, nst=3, nxe, nye 

REAL(iwp):: det, one=1.0_iwp, penalty=1.0e20_iwp, zero—-0.0 iwp 

CHARACTER(LEN=15):: argv, element, dir, type 2d 

}----------- dynamic arrays --———— — ------ 

INTEGER, ALLOCATABLE:: etype(:), g(:), g_g(:, 3: g num(:;,:), kdiag(:), nf(:,:), & 
no(:), node(:), num(:), sense(:) 

REAL(iwp), ALLOCATABLE:: bee(:, :), coord(:, :), dee(:, :), der(:, :), deriv(:, :) & 


eld(:), fun(:), gc(:), g_coord(:, :), jac(:, :), km(:, :), kv(:), loads(:), & 
points(:, UE :), sigma(:), value(:), weights(:), x coords(:), & 
y. coords 


]l----------- input and initialisation -————————-—-— 
CALL getname(argv, nlen) 

OPEN(10, FILE—argv(1:nlen)//".dat") 

OPEN(11, FILE=argv(1:nlen) //’.res’) 

READ(10, *) type 2d, element, nod, dir, nxe, nye, nip, np. types 
CALL mesh size(element, nod, nels, nn, nxe, nye) 
ndof=nod*nodof 


IF(type 2d==‘axisymmetric’) nst=4 
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ALLOCATE(af(nodof, nn), points(nip, ndim), g(ndof), g_coord(ndim, nn), fun(nod), & 
coord(nod, ndim), jac(ndim, ndim), g num(nod, nels), der(ndim, nod), & 
deriv(ndim, nod), bee(nst, ndof), km(ndof, ndof), eld(ndof), weights(nip), & 

g g(ndof, nels), prop(nprops, np types), num(nod), x coords(nxe--1), & 
y. coords(nye--1), etype(nels), gc(ndim), dee(nst, nst), sigma(nst)) 
READ(10, *) prop 

etype=1 

IF(np_types>1) read(10, *) etype 
READ(10, *) x coords, y_coords 

nf=1 
READ(10, *) nr, (k, nf(:, k), i=1, nr) 

CALL formnf(nf) 

neq-MAXVAL(nf) 

ALLOCATE(loads(0:neq), kdiag(neq)) 

kdiag=0 

!--------- loop the elements to find global arrays sizes — — 

elements 1: DO iel=1, nels 
CALL geom rect(element, iel, x. coords, y coords, coord, num, dir) 
CALL num to g(num, nf, g) 

g num(:,iel)2num 

g coord(:, num) TRANSPOSE(coord) 
& g(iel)-g 

CALL fkdiag(kdiag, g) 

END DO elements 1 

CALL mesh(g coord, g num, argv, nlen, 12) 

DO i=2, neq 
kdiag(i)=kdiag(i)+kdiag(i-1) 

END DO 

ALLOCATE(kv(kdiag(neq))) 

WRITE(11, ‘(2(A, I5))) & 

" There are", neq, " equations and the skyline storage is", kdiag(neq) 

l---------- element stiffness integration and assembly — — — 

CALL sample(element, points, weights) 

kv=zero 


gc=one 
elements 2: DO iel=1, nels 

CALL deemat(dee, prop(1, etype(iel)), prop(2, etype(iel))) 

num=g_num(:, iel) 

coord=TRANSPOSE(g_coord(:, num)) 

g=g_g(:, iel) 

km=zero 

int pts 1: DO i=1, nip 
CALL shape fun(fun, points, i) 
CALL shape der(der, points, i) 
jac=MATMUL(der, coord) 
det=determinant (jac) 
CALL invert(jac) 
deriv=MATMUL(jac, der) 
CALL beemat(bee, deriv) 
IF(type_2d==‘axisymmetric’) THEN 

gc=MATMUL(fun, coord) 
bee(4, 1:ndof-1:2)=fun(:) /ge(1) 

END IF 
km=km+MATMUL(MATMUL(TRANSPOSE(bee), dee), bee)*det*weights(i)*ge(1) 

END DO int pts 1 
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CALL fsparv(kv, km, g, kdiag) 
END DO elements_ 2 
loads=zero 
READ(10, *) loaded nodes, (k, loads(nf(:, k)), i=1, loaded nodes) 
READ(10, *) fixed_ freedoms 
IF (fixed_freedoms/=0) THEN 
ALLOCATE(node(fixed_ freedoms), sense(fixed freedoms), & 
value(fixed freedoms), no(fixed freedoms)) 
READ(10, *) (node(i), sense(i), value(i), i=1, fixed_ freedoms) 
DO i=1, fixed freedoms 
no(i)=nf(sense(i), node(i)) 
END DO 
kv(kdiag(no))—kv(kdiag(no))--penalty 
loads(no)=kv(kdiag(no)) *value 
END IF 
l----------- equation solution — - - - - - - - - - - - — 
CALL sparin(kv, kdiag) 
CALL spabac(kv, loads, kdiag) 
loads(0)=zero 
IF(type_2d==‘axisymmetric’) THEN 
WRITE(11, ‘(/A)’)"Noder-dispz-disp" 
ELSE 
WRITE(11, ‘(/A)’)"Nodex-dispy-disp" 
END IF 
DO k=1, nn 
WRITE(11, ‘(I5, 2E12.4)’) k, loads(nf(:, k)) 
END DO 
l----------- recover stresses at nip integrating points — — — 
WRITE(11, ‘(/A, I2, A)’) " The integration point (nip=", nip, ") stresses are:" 
IF(type_2d==‘axisymmetric’) THEN 


WRITE(11, ‘(A, A)’)" Element r-coord z-coord", & 
Y sig_r sig z tau rz sig t" 
ELSE 
WRITE(11, ‘(A, A)’) " Element x-coordy-coord", & 
i sig x sig y tau xy" 
END IF 


elements 3: DO iel=1, nels 
CALL deemat(dee, prop(1, etype(iel)), prop(2, etype(iel))) 
num=g_num(:, iel) 
coord=TRANSPOSE(g_coord(:, num)) 
g=g_g(:, iel) 
eld=loads(g) 
int_pts_2: DO i=1, nip 
CALL shape fun(fun, points, i) 
CALL shape der(der, points, i) 
gc-MATMUL(fun, coord) 
jac=MATMUL(der, coord) 
CALL invert (jac) 
deriv=MATMUL(jac, der) 
CALL beemat(bee, deriv) 
IF(type_2d==‘axisymmetric’) THEN 
gc=MATMUL(fun, coord) 
bee(4, 1:ndof-1:2)=fun(:) /ge(1) 
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END IF 

sigma=MATMUL(dee, MATMUL(bee, eld)) 

WRITE(11, ‘(15, 6E12.4)’) iel, gc, sigma 

END DO int pts 2 

END DO elements 3 
I----------- recover stresses at node points — — — 
points(1,1)—-1, points(1,2)=1, points(2,1)=1, points(2,2)=1 
points(3,1)=-1, points(3,2)=-1, points(4,1)=1, points(4,2)=-1 
WRITE(11,'(/A,12,A)')" The Node point stresses are:" 
IF(type 2d—-'axisymmetric') THEN 


WRITE(11,(A,A))" Element r-coord — z-coord", & 
r sig r sig z w un — qs (UU 

ELSE 
WRITE(11,(A,A))" Element x-coord ^ y-coord", & 
"sig x sig y tau xy" 

END IF 


elements 4: DO iel=1,nels 
CALL deemat(dee,prop(1,etype(iel)) ,prop(2,etype(iel))) 
num=g_num(:,iel) 
coord=TRANSPOSE(g_coord(:,num)) 
g=g_g(:jiel) 
eld=loads(g) 
int pts 3: DO i=1,nip 

CALL shape_fun(fun,points,i) 

CALL shape_der(der,points,i) 

gc=MATMUL(fun,coord) 

jac=MATMUL(der,coord) 

CALL invert (jac) 

deriv=MATMUL(jac,der) 

CALL beemat (bee,deriv) 

IF(type 2d—-'axisymmetric') THEN 
gc-MATMUL(fun,coord) 
bee(4,1:ndof-1:2)=fun(:) /ge(1) 

END IF 

sigma=MATMUL(dee,MATMUL(bee,eld)) 

WRITE(11,'(15,6E12.4)')iel,gc,sigma 

END DO int pts 3 
END DO elements 4 
CALL dismsh(loads, nf, 0.05. iwp, g coord, g num, argv, nlen, 13) 
CALL vecmsh(loads, nf, 0.05. iwp, 0.1 iwp,g coord, g num, argv, nlen, 14) 
STOP 
END PROGRAM p2DFEM 


9.2.2 Numerical Example 


'The input parameters for the above geometric model and the basic conditions are 
contained in file 2DFEM.dat; two elements with uniform size are used; the input 
data and their interpretation are listed in table 9.2. 
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TAB. 9.2 — Input data and interpretation of the two-dimensional program of plane strain 


problem. 


Input data 


‘plane’ 


‘quadrilateral’ 4 ‘y’ 


2141 
1.0e6 0.3 

0.0 15.0 30.0 
0.0 -10.0 

5 


100 200 400 500 600 


1 
3 0 -1.0e6 
0 


Interpretation 

element type 

element, node, direction 
element in z, element in y, integration points, material 
E v 

X coords 

y coords 

nr displacement BC 

k, nf(:,k), i=1, nr 

nr traction BC 

k, nf(:,k), i=1, nr 

fixed freedoms 


Based on running the two-dimensional program of the plane strain problem, the 
output solutions for the above geometric model and basic conditions are contained 
in file 2DFEM.res, and are shown below: 


Node x-disp 
1 0.0000E4 
2 0.0000E- 
3 0.5855E- 
4 0.0000EJ 
5 0.0000E- 
6 0.0000E- 


Element x-coord 
1 0.3170E+01 -0.2113E4 


NNN 


0.1183E4 
0.3170EA 
0.1183EJ 
0.1817 EJ 
0.2683E4 
0.1817E4 
2 0.2683E4 


y-disp 


+00 0.0000E--00 
+00 0.0000E--00 


16 -0.6592E+00 


+00 0.0000E+00 
+00 0.0000E+00 
+00 0.0000E+00 
The integration point (nip= 4) stresses are: 
The integration point (nip= 4) stresses are: 


y-coord 


The node point stresses are: 


Element x-co 

1 0.0000E4 
0.1500E+02 
0.0000E4-00 
0.1500E+02 
0.1500E+02 
0.3000E+02 
0.1500E+02 
0.3000E+02 


NONNNR RR 


ord 
+00 


y-coord 


-0.1000E+ 
-0.1000EA 


-0.1000EA 


-0.1000EA 


sig x sig y tau Xy 


-01 -0.8036E--04 -0.1875E+05 -0.1333E+05 
-02 -0.2113E+01 -0.2999E 4-05 -0.6998E--05 -0.1333E+05 
F01 -0.7887E+01 -0.8036E 4-04 -0.1875E+05 -0.3572E--04 
+02 -0.7887E+01 -0.2999E 4-05 -0.6998E+05 -0.3572E--04 
-02 -0.2113E+01 -0.2999E+05 -0.6998E--05 0.1333E+05 
-02 -0.2113E+01 -0.8036E+04 -0.1875E+05 0.1333E+05 
+02 -0.7887E+01 -0.2999E+05 -0.6998E+05 0.3572E--04 
+02 -0.7887E+01 -0.8036E+04 -0.1875E+05 0.3572E--04 


sig x sig y tau Xy 


0.0000E+00 0.5255E-11 0.2252E-11 -0.1690E--05 

0.0000E--00 -0.3803E 4-05 -0.8873E--05 -0.1690E+05 
-02 0.0000E+00 0.0000E--00 0.0000E+00 
+02 -0.3803E--05 -0.8873E+05 0.2252E-11 
0.0000E--00 -0.3803E+05 -0.8873E+05 0.1690E+05 
0.0000E+00 -0.5255E-11 -0.2252E-11 0.1690E--05 

-02 -0.3803E--05 -0.8873E+05 0.2252bE-11 
-02 0.0000E+00 0.0000E+00 0.0000E+00 
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9.2.3 Interactive Interface 


Figure 9.18 shows that this interface includes the parameter input window on the 
upper left, the model diagram on the upper right, the summarised input data on the 
lower left, and the output data on the lower right. 


& Two-dimensional program of plane strain problem = D x 
Element type 
Element Node Direction 
Element number in x Element number in y Integration points Material 
Young's modulus E Poisson's ratio v | —— — R 
x-coordinate of node Add E F| 
y-coordinate of node | o| 1 3 5 
Dis. BC no. 7 " 
q Q b 
Node number of Dis. BC Value of Dis. BC Add 
Tra. BC no. 
Node number of Tra. BC Value of Tra. BC Add a 
Fixed freedom 
Display input data Generate file 2DFEM. dat Compute Element number: © and © Node number: 1-6 
Dry cupa dum ppm B Model of two-dimensional program of plane strain problem 
Display input data: Display output data: 
D—————————————————————————— 


Fic. 9.18 — Interface of the two-dimensional program of plane strain problem. 


(1) Step 1 


Figure 9.19 shows one parameter in the first line of the parameter input window: 
Element type: element type of two-dimensional elements used. The value entered 
in this figure was ‘plane’. 


(2) Step 2 


Figure 9.20 shows three parameters in the second line of the parameter input 
window: 

Element: element category in the selected element type. The value entered in this 
figure was ‘quadrilateral’. 

Node: node number of the element. The value entered in this figure was 4. 

Direction: direction in which elements are numbered. The value entered in this 
figure was ‘y’. 


(3) Step 3 


Figure 9.21 shows four parameters in the third line of the parameter input 
window: 
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(Two-dimensional program of plane strain problem - o x 


Node Direction 
Element number in y Integration points Material 


Node number of Dis. BC Value of Dis. BC od 


Node number of Tra. BC Value of Tra. BC Add 


Display input data Generate file 2DFEM dat Compute Element number: © and © Node number: 1-6 
" z Model of two-dimensional program of plane strain problem 
Display output data Retum Exit 


Display input data: Display output data: 


Fic. 9.19 — Element type in the two-dimensional program of plane strain problem. 


f. Two-dimensional program of plane strain problem - 0 x 
Element type [plane — 
Element fauadrilat Node + Direction 4 
Element number in x Element number in y Integration points Material 
Young's modulus E Possonsratiov | —— — 
x-coordinate ofnode | _Add_| 
y-coordinate ofnode | — Bd 
Dis. BC no. [| 
Node number of Dis. BC =E Value of Dis. BC Add 
Tra. BC no. mum 
Node number of Tra. BC [ Value of Tra. BC Add 
Fixed freedom 
Display input data Generate file 2DFEM dat Compute Element number: @ and @) Node number: 1-6 
Dipl ompar dui pum Exit Model of two-dimensional program of plane strain problem 


[Display input data: 


Fic. 9.20 — Element category, node, and direction in two-dimensional program of plane strain 
problem. 
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I5. Two-dimensional program of plane strain problem - n x 

pes plane 

Element quadrilat — Node 4 Direction 2 

Element number in x |2 Element number in y |1 Integration points |4 

Young's modulus E Poisson's ratio v 

x-coordinate of node a 

y-coordinate of node Add 

Dis. BC no. 
Node number of Dis. BC Value of Dis. BC Add $ 
Tra. BC no. 
Node mmber of Tra. BC Value of Tra BC Add = 
Fixed freedom b a € 

Display input data Generate fle 2DFEM dat Compute Element number: © and © Node number: 1-6 
=e a = Model of two-dimensional program of plane strain problem 

Display input data: Display output data: 


Fic. 9.21 — Element number in z, y, integration points, and material in the two-dimensional 
program of plane strain problem. 


Element number in a: element number in z-direction. The value entered in this 
figure was 2. 

Element number in y: element number in y-direction. The value entered in this 
figure was 1. 

Integration points: number of integration points. The value entered in this figure 
was 4. 

Material: number of material types. The value entered in this figure was 1. 


(4) Step 4 


Figure 9.22 shows two parameters in the fourth line of the parameter input 
window: 

Young's modulus E: Young's modulus. The value entered in this figure was 
1.0e6. 

Poisson's ratio v: Poisson's ratio. The value entered in this figure was 0.3. 


(5) Step 5 


Figure 9.23 shows one parameter in the fifth line of the parameter input window: 

az-coordinate of node: the z-coordinate of each node. The z-coordinate value of 
each node can be input sequentially, and the option ‘Add’ is selected. Subsequently, 
the window *Add successfully' appears, and the option *Agree' is selected to confirm 
the addition. The values entered in this figure were 0, 15, and 30. 
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P3. Two-dimensional program of plane strain problem 
Element type [plane 

Element [auadrilat Node 

Element number in x |2 Element number in y 
Young's modulus E. 1.0e6 Poisson's ratio v 
x-coordinate of node 
y-coordinate of node 
Dis. BC no. 


Node number of Dis. BC 
Tra. BC no. 


Node number of Tra. BC Value of Tra. BC 


Fixed freedom 


Display input data. Generate file 2DFEM dat Element number: @ and © Node number: 1-6 
* r Model of two-dimensional program of plane strain problem 
Display output data Retum 


[Display input data: Display output data: 


Fic. 9.22 — Young’s modulus and Poisson’s ratio in the two-dimensional program of plane 
strain problem. 


plane — 
[quadrilat 
Element number in x — 2 
Young's modulus E 
x-coordinate of node 


y-coordinate of node 
Dis. BC no. 

Node number of Dis. BC 
Tra. BC no. 

Node number of Tra. BC 
Fixed freedom 


Display input data en Element number: © and © Node number: 1-6 


Model of two-dimensional program of plane strain problem 


Display output data 


Display input data: Display output data: 


Fic. 9.23 — The xcoordinate of node in the two-dimensional program of plane strain 
problem. 


138 Basic Theory of Finite Element Method 


(6) Step 6 


Figure 9.24 shows one parameter in the sixth line of the parameter input window: 

y-coordinate of node: The ycoordinate of each node. The ycoordinate value of 
each node can be input sequentially, and the option ‘Add’ is selected. Then, the 
window ‘Add successfully’ appears, and the option ‘Agree’ is selected to confirm the 
addition. The values entered in this figure were 0 and —10. 


5 


Element type plane 


e quadrilat Node 4 Direction y 
Element numberin x |2 Element number in y | Integration points |4 
Young's modulus E 1.0e6 Poisson's ratio v [o3 ——— 
x-coordinate ofnode — |30 . Adi | 
y-coordinate ofnode |0 
Dis. BC no. 
Node mumber of Dis. BC Value of Dis. BC Add | d 
Tra. BC no. — 
Node number of Tra. BC yag owe — 0 X us E 
Fixed freedom 
Display input data 444 vocctestuly CES Element number: © and © Node number: 1-6 
eee! === E Model of two-dimensional program of plane strain problem 


Display input data: Display output data: 


Fic. 9.24 — The y-coordinate of node in two-dimensional program of plane strain problem. 


(7) Step 7 


Figure 9.25 shows one parameter in the seventh line of the parameter input 
window: 

Dis. BC no.: the number of displacement boundary conditions. The value entered 
in this figure was 5. 


(8) Step 8 


Figure 9.26 shows two parameters in the eighth line of the parameter input 
window: 

Node number of Dis. BC: node number of displacement boundary, specifying the 
node at which a specific displacement occurs. The values entered in this figure were 
1, 2, 4, 5, and 6. 
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5. Two-dimensional program of plene strain problem - o0 x 
Element type [plane — 
Element jquadrilat Node 4 Direction b 
Element number in x — |? Element number in y |t Integration points [# —— 
Young's modulus E 1.0e6 Poisson's ratio v [9.3 

x-coordinate ofnode |30 __ Add] 

y-coordinate ofnode — [10 LESS] 
Dis. BC no. 5 
Node number of Dis. BC Value of Dis.BC | —— Add P 
Tra. BC no. 
Node mmber of Tra. BC Value of Tra BC Add ~ 
Fixed freedom 

Display input data Generate file 2DFEM dat Compute Element number: © and © Node number: 1-6 

DUpisy oupu dus Prem Exit Model of two-dimensional program of plane strain problem 

[Displayinputdatas SS Display output data: 


Fic. 9.25 — Number of displacement boundary conditions in the two-dimensional program of 
plane strain problem. 


Bi p — B 

Element type plane 
Element quadrila! — Node 4 Direction y 
Element numberin x 2 Element number in y |! Integration points |4 Material |1 
Young's modulus E 1.0e6 Poisson's ratio 

'oung's 'oisson's ratio v [0.3 P 
x-coordinate of node 30 Add x F | 
y-coordinate of node -10 Add oj 1 3 5 
Dis. BC no. s di " 
Node number of Dis. BC |1 Value of Dis. BC |0 o Add 2, 
Tra. BC no. 

— STwedme. — O X B 

Node number of Tra. BC Value of Tra. L1 
Fixed freedom 

Display input data Generate file te | Element number: © and © Node number: 1-6 

Disp dus LI Model of two-dimensional program of plane strain problem 


[Display input data: Display output data: 


Fic. 9.26 — Node number and value of displacement boundary condition in the 
two-dimensional program of plane strain problem. 
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Value of Dis. BC: values of displacement boundary, specifying the values used for 
displacement boundary in z- and y-directions, respectively. The value of the dis- 
placement boundary of each specifying node can be input sequentially, and click the 
option ‘Add’ is selected. Then, the window ‘Add successfully’ appears, and the 
option ‘Agree’ is selected to confirm the addition. The values entered in this figure 
were 0 and 0. 


(9) Step 9 


Figure 9.27 shows one parameter in the ninth line of the parameter input 
window: 

Tra. BC no.: the number of traction boundary conditions. The value entered in 
this figure was 1. 


© Two-dimensional program of plane strain problem - n x 

Element type plane 

Element iquadrilat Node 4 Direction y 

Element number in x — |? Element number in y |! Integration points |4 Material |1 

Young's modulus E 1.0e6 Poisson's ratio v [0.3 

x-coordinate of node [30 ae 

y-coordinate ofnode |-10 ae | 
Dis. BC no. 5 
Node mmber of Dis. BC |6 Value of Dis. BC |0 o LES] ? 
Tra. BC no. - 
Node number of Tra. BC Value of Tra BC Add B 
Fixed freedom 

Display input data Generate file 2DFEM dat Compute Element number: © and © Node number: 1-6 

=a hes Ext Model of two-dimensional program of plane strain problem 
Display input data: Display output data: 


Fic. 9.27 — Number of traction boundary conditions in the two-dimensional program of plane 
strain problem. 


(10) Step 10 


Figure 9.28 shows two parameters in the tenth line of the parameter input 
window: 
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B 


Element type plane 
Element quadrilat ^ Node 4 
Element numberin x |2 Element number in y |1 
Young's modulus E 1.0e6 Poisson's ratio v [o3 — 
x-coordinate of node 3o Add 
y-coordinate of node -10 Add 
Dis. BC no. ü 
Node number of Dis. BC |6 Value of Dis. BC e 
Tra. BC no. 1 
Node number of Tra. BC |3 Value of Tra BC ~ 
Fixed freedom 
Display input data Generate file 2DÉ Element number: @ and © Node number: 1-6 
meaa] ma Model of two-dimensional program of plane strain problem 


Display input data: isplay output data: 


Fic. 9.28 — Node number and value of traction boundary condition in the two-dimensional 
program of plane strain problem. 


Node number of Tra. BC: node number of traction boundary, specifying the 
node at which specific traction occurs. The value entered in this figure was 3. 

Value of Tra. BC: values of traction boundary, specifying the values used for the 
traction boundary in z- and ydirections, respectively. The value of traction 
boundary of each specifying node can be input sequentially, and the option 
‘Add’ is selected. Then, the window ‘Add successfully’ appears, and the option 
* Agree' is selected to confirm the addition. The value entered in this figure were 0 
and —1.0e6. 


(11) Step 11 


Figure 9.29 shows one parameter in the eleventh line of the parameter input 
window: 

Fixed freedom: the number of nodes with the specified displacement. The value 
entered in this figure was 0 by default. 


(12) Step 12 


As shown in figure 9.30, the option ‘Display input data’ is selected, and the 
aforementioned input dataset is summarised on the lower left of the interface. 
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© Two-dimensional program of plane strain problem 


Element type plane 
fice [quadrilat 


E 
UN 
2m 
= 


Element number in x 
Young's modulus E. 
x-coordinate of node 


y-coordinate of node 

Dis. BC no. 5 
Node number of Dis. BC |6 
Tra. BC no. 1 


Node number of Tra. BC |3 


Node j4 
Element number in y |1 
Poisson's ratio v [ns —— 


Add 


Add 


Value of Dis. BC |0 
Value of Tra. BC |0 


[3 


-1.0e6 


Element number: © and © Node number: 1~6 


Model of two-dimensional program of plane strain problem 


[Display input data: 


Fic. 9.29 — Fixed freedom 


Display output data: 


in the two-dimensional program of plane strain problem. 


I Two-dimensional program of plane strain problem 


Element type plane 
Element 

Element number in x 
Young's modulus E. 
x-coordinate of node 


iquadrilat 


y-coordinate of node 
Dis. BC no. 

Node number of Dis. BC 
Tra. BC no. 

Node number of Tra. BC 


Fixed freedom 


(7 ] 
Display output data 


Node 
Element number in y 


Poisson's ratio v 


Value of Tra. BC 


Rem | 


Element number: @ and (2) Node number: 1~6 


Model of two-dimensional program of plane strain problem 


plane" 
quadrilateral’ 4 'y' 
2141 
1.0e6 0.3 

15 30 

-10 
5 
100200400500600 
1 
3 0 -1000000 


Display output data: 


Fic. 9.30 — Input data display in the two-dimensional program of plane strain problem. 
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(13) Step 13 

As shown in figure 9.31, the option ‘Generate file 2DFEM.dat’ is selected, and 
the input dataset is recorded in the file 2DFEM.dat, which is saved in the path of the 
current file. Subsequently, the window ‘Generate successfully’ appears, and the 
option ‘Agree’ is selected for confirmation. 


5 
Element type plane 
Element quadrilat Node a 
Element number in x — |2 feba un». 
Young's moduhis E  [L0e6 Poisson's ratio v [o3 — 
x-coordinate ofnode |30 Add 
y-coordinate ofnode — 10 Add 
Dis. BC no. 5 
Node number of Dis. BC |6 Value of Dis. BC |0 b 
Tra. BC no. 1 
Node number of Tra. BC |3 Value of Tra. BC) © Two-dime... = 
Fixed freedom o 
Display put cian tl Element number: © and © Node number: 1-6 
Display output data Ll Model of two-dimensional program of plane strain problem 


aoe ' isplay output data: 
[quadrilateral' 4 'y' 
2141 
1.0e6 0.3 

15 30 

-10 


I5 
100200400500600 


" 
3 0 -1000000 
jo 


Fic. 9.31 — Input data generation in the two-dimensional program of plane strain problem. 


(14) Step 14 


As shown in figure 9.32, the option ‘Compute’ is selected, and the computation 
program starts to conduct operations using the aforementioned input data. Once 
the computation is completed, the window *Compute successfully? appears, and the 
option *Agree' is selected for confirmation. The output dataset above is recorded in 
the file 2DFEM.res, which was saved in the path of the current file. 


(15) Step 15 


As shown in figure 9.33, the option ‘Display output data’ is selected, and the 
output data of computed results are displayed on the lower right of the interface. 

The option ‘Return’ in the interactive interface, will return to the interface for 
selecting the computation program shown in figure 9.2. The option 'Exit will exit 
the current interactive interface. 
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plane 


Element type 
Element 

Element number in x 
Young's modulus E. 
x-coordinate of node 
y-coordinate of node 
Dis. BC no. 


quadrilat 


Node number of Dis. BC 
Tra. BC no. 

Node number of Tra. BC 
Fixed freedom. 


Display input data 
Display output data 


Value of Dis. BC 


Value of Tra. BQ 


Generate file 2D 
ELE 


B Two-dime.. — [s] 


Compute successfully 


o [o __Add 


x 


Element number: © and © Node number: 1-6 


Model of two-dimensional program of plane strain problem 


‘plane’ 
‘quadrilateral’ 4 'y' 
2141 

1.0e6 0.3 

[0 15 30 

jo -10 

I 
100200400500600 
" 
|3 0 -1000000 
jo 


Display output data: 


Fic. 9.32 — Computation control in the two-dimensional program of plane strain problem. 


f. Two-dimensional program of plane strain problem 


Element type 

Element 

Element number in x 
Young's modulus E 
x-coordinate of node 
y-coordinate of node 
Dis. BC no. 

Node number of Dis. BC 
Tra. BC no. 

Node number of Tra. BC 


Fixed freedom. 
Display input data 
[ put dan | 


[e| 
=] 


a 


Element number: © and (2) Node number: 1-6 


Model of two-dimensional program of plane strain problem 


[plane 
‘quadrilateral’ 4 'y’ 
2141 

1.0e6 0.3 

jo 15 30 

jo -10 

5 
100200400500600 
" 
[3 0 -1000000 
jo 


Node xdisp  y-disp 

1 0.0000E+00 0.0000E+00 
2 0.0000E+00 0.0000E+00 
3 0.5855E-16 -0.6592E+00 
4 0.0000E+00 0.0000E+00 
5 0.0000E+00 0.0000E+00 
6 0.0000E+00 0.0000E+00 


The integration point (nip= 4) stresses are: 

Elementx-coord y-coord sigx sigy  tauxy 
1 0.3170E+01 -0.2113E+01 -0.8036E+04 -0.1875E«05 -0.1333E+05 
1 0.1183E+02 -0.2113E+01 -0.2999E+05 -0.6998E+05 -0.1333E+05 


Fic. 9.33 — Output data display in the two-dimensional program of plane strain problem. 
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9.3 Three-Dimensional Program of Solid Compression 


'The geometric model of three-dimensional program of solid compression is shown in 
figure 9.34, with dimensions a = 0.5 in the horizontal z-direction, b = 10 in the 
horizontal y-direction, and c = 2 in the vertical z-direction. The material parameters 
are Young's modulus E = 1 X 10° and Poisson's ratio v = 0.3. The boundary 
conditions are that the upper left and upper right vertices are fixed, and the lower 
parts ao fixed. The vertical concentrated load F acting on the middle domain is 
-1 x 10°. 


Element number: (Dand() 
Node number: 1~12 


Fic. 9.34 — Geometric model of the three-dimensional program of solid compression. 


9.3.1 Main Program 


The main program of the three-dimensional program of solid compression is shown 
below and contains the following modules: dynamic arrays, input, and initialization, 
loop the elements to find global arrays sizes, element stiffness integration and 
assembly, equation solution, recover stresses at nip integrating points and recover 
stresses at node points. 
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PROGRAM p3DFEM 

USE main 

USE geom 

IMPLICIT NONE 

INTEGER, PARAMETER:: iwp=SELECTED_REAL_KIND(15) 

INTEGER:: fixed freedoms, i, iel, k, loaded nodes, ndim=3, ndof, nels, neq, nip, & 
nlen, nn, nprops=2, np types, nod, nodof=3, nr, nst—6, nxe, nye, nze 

REAL(iwp):: det, penalty=1.0e20_iwp, zero=0.0_iwp 

CHARACTER(LEN=15):: argv, element=‘hexahedron’ 

!----------- dynamic arrays -————-—--------- 

INTEGER, ALLOCATABLE:: etype(:), g(:), g_g(:, :), g_num(:, :), kdiag(:), nf(:, :), & 
no(:), node(:), num(:), sense(:) 

REAL(iwp), ALLOCATABLE:: bee(:, :), coord(:, :), dee(:, :), der(:, :), deriv(:,:), & 
eld(:), fun(:), gc(:), g_coord(:, :), jac(:, :), km(:, :), kv(:), loads(:), & 
points(:, :), prop(:, :), sigma(:), value(:), weights(:), x coords(:), & 
y_coords(:), z coords(:) 

!----------- input and initialisation - ———————--— 

CALL getname(argv, nlen) 

OPEN(10, FILE=argv(1:nlen) //’.dat’) 

OPEN(11, FILE=argv(1:nlen) //’.res’) 

READ(10, *) nod, nxe, nye, nze, nip, np types 

CALL mesh size(element, nod, nels, nn, nxe, nye, nze) 

ndof=nod*nodof 

ALLOCATE(nf(nodof, nn), points(nip, ndim), dee(nst, nst), coord(nod, ndim), & 
jac(ndim, ndim), der(ndim, nod), deriv(ndim, nod), g(ndof), bee(nst, ndof), & 
km(ndof, ndof), eld(ndof), sigma(nst), g_g(ndof, nels), g_coord(ndim, nn), & 

g num(nod, nels), weights(nip), num(nod), prop(nprops, np types), & 
x coords(nxe--l), y_coords(nye+1), z_coords(nze+1), etype(nels), fun(nod), & 
gc(ndim)) 

READ(10, *) prop 

etype=1 

IF(np types-1) READ(10, *) etype 

READ(10, *) x coords, y coords,z coords 

nf=1 

READ(10, *) nr, (k, nf(:, k), i=1, nr) 

CALL formnf(nf) 

neq-MAXVAL(nf) 

ALLOCATE(loads(0:neq), kdiag(neq)) 

kdiag=0 

!----------- loop the elements to find global arrays sizes — — — 

elements 1: DO iel=1, nels 
CALL hexahedron xz(iel,x coords, y coords,z coords, coord, num) 
CALL num to g(num, nf, g) 

g num(:,iel)2num 

g coord(:, num) TRANSPOSE(coord) 
g_g(:, iel)=g 

CALL fkdiag(kdiag, g) 

END DO elements_ 1 

DO i=2, neq 
kdiag(i)=kdiag(i)+kdiag(i-1) 

END DO 
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ALLOCATE(kv(kdiag(neq))) 
WRITE(11, ‘(2(A, I5))) & 
" There are", neq, " equations and the skyline storage is", kdiag(neq) 
| —————— —— MÀ element stiffness integration and assembly — — — 
CALL sample(element, points, weights) 
kv=zero 
elements _2: DO iel=1, nels 
CALL deemat(dee, prop(1, etype(iel)), prop(2, etype(iel))) 
num=g_num(:, iel) 
g=g_g(;, iel) 
coord- TRANSPOSE(g coord(:, num)) 
km=zero 
gauss pts 1: DO i=1, nip 
CALL shape_der(der, points, i) 
jac=MATMUL(der, coord) 
det=determinant (jac) 
CALL invert (jac) 
deriv=MATMUL(jac, der) 
CALL beemat(bee, deriv) 
km=km+MATMUL(MATMUL(TRANSPOSE(bee), dee), bee)*det*weights(i) 
END DO gauss pts 1 
CALL fsparv(kv, km, g, kdiag) 
END DO elements 2 
loads=zero 
READ(10, *)loaded nodes 
IF(loaded nodes/—0) READ(10, *) (k, loads(nf(:, k)), i21, loaded nodes) 
READ(10, *) fixed. freedoms 
TF (fixed freedoms/—0) THEN 
ALLOCATE(node(fixed freedoms), sense(fixed freedoms), & 
value(fixed freedoms), no(fixed freedoms)) 
READ(10, *) (node(i), sense(i), value(i), i=1, fixed _ freedoms) 
DO i=1, fixed freedoms 
no(i)=nf(sense(i), node(i)) 
END DO 
kv(kdiag(no))=kv(kdiag(no))+-penalty 
loads(no)=kv(kdiag(no))*value 
END IF 
eee ems equation solution - -———-—---------- 
CALL sparin(kv, kdiag) 
CALL spabac(kv, loads, kdiag) 
loads(0)—zero 
WRITE(11, ‘(/A)’)" Node x-disp y-disp — z-disp" 
DO k=1, nn 
WRITE(11, *(I5, 3E12.4)’) k, loads(nf(:, k)) 
END DO 


recover stresses at nip integrating points — — — 
WRITE(11, ‘(/A, I2, A)’)" The integration point (nip=", nip, ") stresses are:" 
WRITE(11, ‘(A, /, A)’)" | Element x-coord y-coord z-coord" & 
"sig x sig y sig z tau xy tau yz tau zx" 
elements 3: DO iel=1, nels 

CALL deemat(dee, prop(1, etype(iel)), prop(2, etype(iel))) 

num=g_num(:, iel) 

coord=TRANSPOSE(g_coord(:, num)) 

g=g_g(:, iel) 
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eld=loads(g) 
gauss pts 2: DO i=1, nip 
CALL shape_der(der, points, i) 
CALL shape fun(fun, points, i) 
gc=MATMUL(fun, coord) 
jac=MATMUL(der, coord) 
CALL invert(jac) 
deriv MATMUL(jac, der) 
CALL beemat(bee, deriv) 
sigma=MATMUL(dee, MATMUL(bee, eld)) 
WRITE(11, (I8, 4X, 3E12.4)’) iel, gc 
WRITE(11, ‘(6E12.4)’) sigma 
END DO gauss pts 2 
END DO elements 3 
I----------- recover stresses at node points — — — 


points(1,1)—1, points(2,1)—1, points(3,1)—1, points(4, Md 
points(5,1)=-1, points(6,1)=-1, points(7,1)—-1, points(8,1)— 
points(1,2)—1, points(2,2)—1, Poe 2)=1, pointe " 
points(5,2)=1, points(6,2)=-1, points(7,2)—1, points(8,2)—- 
points(1,3)=1, points(2,3)=-1, points(3,3)—1, points(4,3)—- 
points(5,3)—1, points(6,3)—1, points(7,3) —-1, points(8,3)— E 
WRITE(11,'(/A,12,A)')" The Node point stresses are:" 


WRITE(11,(A,/,A))" | Element x-coord y-coord  z-coord" & 
r sig_x sig y sig z tau xy tau yz tau zx" 
elements 4: DO iel=1,nels 
CALL deemat(dee,prop(1,etype(iel)) ,prop(2,etype(iel))) 
num=g_num(:,iel) 
coord=TRANSPOSE(g_coord(:,num)) 
&—g g(iel) 
eld=loads(g) 
gauss pts 3: DO i=1,nip 
CALL shape der(der,points,i) 
CALL shape fun(fun,pointsj) 
c=MATMUL(fun,coord) 
jac=MATMUL(der,coord) 
CALL invert (jac) 
deriv=MATMUL(jac,der) 
CALL beemat(bee,deriv) 
sigma=MATMUL(dee,MATMUL(bee,eld)) 
WRITE(11,(18,4X,3E12.4)))iel,gc 
WRITE(11,'(6E12.4)')sigma 
END DO gauss pts 3 
END DO elements 4 
STOP 
END PROGRAM p3DFEM 


9.3.2 Numerical Example 


The input parameters for the above geometric model and basic conditions are 
contained in file 3DFEM.dat; two elements with uniform size are used; the input 
data and their interpretation are listed in table 9.3. 
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TAB. 9.8 — Input data and interpretation of the three-dimensional program of solid 


compression. 

Input data Interpretation 

8 element-node 

12181 element in z, element in y, element in z, 
integration points, material 

1.0e6 0.3 E v 

0.0 0.5 € coords 

0.0 1.5 3.0 y. coords 

0.0 -2.0 z_coords 

10 nr displacement BC 

1000 2000 3000 4000 7000 kon SEL us 

8000 9000 10000 11000 12000 

2 nr traction BC 

5 0.0 0.0 -1.0e6 6 0.0 0.0 -1.0e6 k, nf(:,k), i=1, nr 

0 fixed_ freedoms 


Based on running the three-dimensional program of solid compression, the 
output solutions for the above geometric model and basic conditions are contained 
in file 3DFEM.res, and are shown below: 


Node x-disp 
0.0000E+00 0.0000E- 
0.0000E+00 0.0000E- 
0.0000E+00 0.0000E- 


= 


y-disp 


z-disp 
-00 0.0000E- 
-00 0.0000E- 
-00 0.0000E- 


-00 
-00 
-00 


0.0000E+00 0.0000E+00 0.0000E+00 
-0.3438E--00 0.1201E-15 -0.4332E+01 
0.3438E4-00 0.1728E-15 -0.4332E+01 
0.0000E+00 0.0000E+00 0.0000E-+00 
0.0000E+00 0.0000E+00 0.0000E-+00 
9 0.0000E+00 0.0000E+00 0.0000E+00 
10 0.0000E+00 0.0000E--00 0.0000E4-00 
11 0.0000E4-00 0.0000E--00 0.0000E4-00 
12 0.0000E+00 0.0000E--00 0.0000E4-00 
The integration point (nip— 8) stresses are: 


oo -10» Ct à C2 F2 


Element | x-coord y-coord z-coord 
sig x sig y sig Z tau Xy tau yz tau zx 
1 0.3943E--00 0.1183E+01 -0.4226E+00 


0.1661E+02 -0.4920E+02-0.1806E+03 0.4014E+01-0.8760E+02 0.3011E+01 
1 0.3943E--00 0.1183E+01 -0.1577E+01 

-0.6770E+02-0.8533E+02-0.2167E+03 0.1076E--01-0.2347E--02 0.3011E+01 
1 0.3943E+00 0.3170E+00 -0.4226E+00 

0.4450E--01 -0.1318E+02-0.4839E+02 0.4014E+01 -0.87600E+02 0.8068E+00 
1 0.3943E+00 0.3170E+00 -0.1577E+01 

-0.1814E+02-0.2286E+02-0.5807E+02 0.1076E+01-0.2347E+02 0.8068E--00 
1 0.1057E--00 0.1183E+01 -0.4226E+00 
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0.1661E--02 -0.4920E+02 -0.1806E+03 -0.4014E+01 -0.8760E 4-02 -0.3011E+01 


1 0.1057E--00 0.3170E+00 -0.4226E-4-00 


0.4450E--01 -0.1318E4-02 -0.4839E4-02 -0.4014E4-01 -0.8760E 4-02 -0.8068E+00 


T 0.1057E--00 0.1183E+01 -0.1577E+01 


-0.6770E--02 -0.8533E+02 -0.2167E+03 -0.1076E+01 -0.2347E+02 -0.3011E+01 


1 0.1057E--00 0.3170E+00 -0.1577E+01 


-0.1814E+02 -0.2286E+02 -0.5807E+02 -0.1076E+01 -0.2347E+02 -0.8068E+00 


2 0.3943E--00 0.2683E+01 -0.4226E+00 


0.4450E--01 -0.1318E+02 -0.4839E+02-0.4014E+01 0.8760E+02 0.8068E+00 


2 0.3943E--00 0.2683E+01 -0.1577E+01 


-0.1814E+02-0.2286E+02-0.5807E+02-0.1076E+01 0.2347E+02 0.8068E--00 


2 0.3943E--00 0.1817E+01 -0.4226E+00 


0.1661E--02 -0.4920E+02-0.1806E+03-0.4014E+01 0.8760E+02 0.3011E--01 


2 0.3943E--00 0.1817E+01 -0.1577E+01 


-0.6770E+02-0.8533E+02-0.2167E+03-0.1076E+01 0.2347E+02 0.3011E+01 


2 0.1057E--00 0.2683E+01 -0.4226E+00 


0.4450E--01-0.1318E4-02-0.4839E4-02 0.4014E+01 0.8760E+02-0.8068E+00 


2 0.1057E+00 0.1817E+01 -0.4226E-+00 


0.1661E--02-0.4920E--02-0.1806E--03 0.4014E+01 0.8760E+02-0.3011E+01 


2 0.1057E--00 0.2683E+01 -0.1577E+01 


-0.1814E+02-0.2286E+02-0.5807E+02 0.1076E--01 0.2347E+02-0.8068E+00 


2 0.1057E+00 0.1817E+01 -0.1577E+01 


-0.6770E+02-0.8533E+02-0.2167E+03 0.1076E--01 0.2347E+02-0.3011E+01 


The node point stresses are: 
Element x-coord y-coord z-coord 


sig x sig y sig Z tau Xy tau yz tau zx 


1  0.5000E4-00 0.1500E+01 0.0000E--00 


0.6018E+02 -0.4561E+02-0.2122E+03 0.8816E--01-0.1111E--03 0.6612E+01 


1  0.5000E4-00 0.1500E--01 -0.2000E+01 


-0.1250E--03 -0.1250E+03 -0.2916E+03 0.0000E+00 0.3323E-14 0.6612E+01 


1 -0.5000E+00 0.0000E+00 0.0000E--00 

0.6645E-14 0.1551E-13 0.6645E-14 0.8816E+01 -0.1111E--03 0.0000E--0 
1  0.5000E-4-00 0.0000E--00 -0.2000E+01 

0.0000E--00 0.0000E+00 0.0000E4-00 0.0000E+00 0.0000E--00 0.0000E4 
1  0.0000E4-00 0.1500E+01 0.0000E--00 


0 


-00 


0.6018E--02 -0.4561E 4-02 -0.2122E+03 -0.8816E--01 -0.1111E 4-03 -0.6612E+01 


1  0.0000E4-00 0.0000E+00 0.0000E--00 


0.4618E-14 0.1078E-13 0.4618E-14 -0.8816E+01 -0.1111E--03 0.0000E+00 


1  0.0000E-4-00 0.1500E--01 -0.2000E+01 


-0.1250E--03 -0.1250E+03 -0.2916E--03 0.0000E--00 0.2309E-14 -0.6612E+01 


1  0.0000E4-00 0.0000E--00 -0.2000E+01 
0.0000E--00 0.0000E+00 0.0000E4-00 0.0000E4-00 0.0000E--00 0.0000E4 
2  0.5000E+00 0.3000E+01 0.0000E4-00 


-00 


-0.6645E-14 -0.1551E-13 -0.6645E-14 -0.8816E+01 0.1111E--03 0.0000E4-00 


2  0.5000E--00 0.3000E+01 -0.2000E4-01 


0.0000E--00 0.0000E+00 0.0000E+00 0.0000E--00 0.0000E--00 0.0000E+00 
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-01 0.0000E4-00 


0.6018E--02 -0.4561E+02-0.2122E+03-0.8816E+01 0.1111E+03 0.6612E+01 


2  0.5000E+00 0.1500E- 


-01 -0.2000E--01 


-0.1250E+03 -0.1250E+03 -0.2916E+03 0.0000E+00 0.3323E-14 0.6612E+01 


2 0.0000E+00 0.3000EJ 


-01 0.0000E- 


-00 


-0.4618E-14 -0.1078E-13 -0.4618E-14 0.8816E+01 0.1111E--03 0.0000E+00 


2 0.0000E+00 0.1500E- 


-01 0.0000E- 


+00 


0.6018E--02-0.4561E--02-0.2122E--03 0.8816E--01 0.1111E+03-0.6612E+01 


2  0.0000E+00 0.3000EJ 


+01 -0.2000E--01 


0.0000E--00 0.0000E+00 0.0000E+00 0.0000E--00 0.0000E--00 0.0000E+00 


2 0.0000E4-00 0.1500EJ 


+01 -0.2000E--01 


-0.1250E--03 -0.1250E+03 -0.2916E--03 0.0000E--00 0.2309E-14 -0.6612E+01 


9.39.3 Interactive Interface 


Figure 9.35 shows that this interface includes the parameter input window on the 
upper left, the model diagram on the upper right, the summarised input data on the 


lower left, and the output data on the lower right. 


F3. Three-dimensional program of solid compression 
Element node 
Element number in x 


Young's modulus E Poisson's ratio v 
x-coordinate of node Add 
y-coordinate of node Add 
z-coordinate of node Sad 
Dis. BC no. 
Node number of Dis. BC Value of Dis. BC 
Tra. BC no. 
Node number of Tra. BC Value of Tra. BC 
Fixed freedom 
Display input data Generate file 3DFEM dat 
Display output data Retum 


Element number in y Element number in = Integration points Material 


4 Element number: © and © 
M Node number: 1-12 


Model of three-dimensional program of solid compression 


Display input data: 


Display output data: 


Fic. 9.35 — Interface of the three-dimensional program of solid compression. 


(1) Step 1 


Figure 9.36 shows one parameter in the first line of the parameter input window: 
Element node: node number of the element. The value entered in this figure was 8. 
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E Three-dimensional program of solid compression = Ud Xx 
Element node e 
Element number in x Element number in y Element number in = Integration points Material 
Young's modulus E Poisson's ratio v “Ag 
x-coordinate of node Add a 
y-coordinate of node Add 
z-coordinate of node Add 
Dis. BC no. 
Node number of Dis. BC Value of Dis. BC Add 
Tra. BC no. 
Node number of Tra. BC Value of Tra. BC ES $ 2N 
Fixed freedom 

Display input data Generate file 3DFEM dat Compute 
Element number: (1) and (2) 
Display output data Return Exit ls Node number: 1-12 

IEEE -—( Model of three-dimensional program of solid compression 
[Display input data: Display output data: 


Fic. 9.36 — Element node in the three-dimensional program of solid compression. 


(2) Step 2 


Figure 9.37 shows five parameters in the second line of the parameter input 
window: 

Element number in a: element number in z-direction. The value entered in this 
figure was 1. 

Element number in y: element number in y-direction. The value entered in this 
figure was 2. 

Element number in z: element number in zdirection. The value entered in this 
figure was 1. 

Integration points: number of integration points. The value entered in this figure 
was 8. 

Material: number of material types. The value entered in this figure was 1. 


(3) Step 3 


Figure 9.38 shows two parameters in the third line of the parameter input 
window: 

Young's modulus E: Young's modulus. The value entered in this figure was 
1.0e6. 

Poisson's ratio v: Poisson's ratio. The value entered in this figure was 0.3. 
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© Three-dimensional program of solid compression - o X 


Node number of Tra. BC 
Fixed freedom. 


Element number: © and © 
Node number: 1-12 


Model of three-dimensional program of solid compression 


Display input data: Display output data: 


Fic. 9.37 — Element number in z, y, z, integration points, and material in the three- 
dimensional program of solid compression. 


B Three-dimensional program of solid compression = o x 
Element node 8 
Element number in x 1 Element number in y |? Element number in = Integration points |8 Material |1 
Vomgemodske K, 1.0e6 Béksoneraio v. 03 As 
x-coordinate of node Add e 
y-coordinate of node Add 
z-coordinate of node oa 
Dis. BC no. Em 
Node number of Dis. BC Value of Dis. BC oe o 
Tra. BC no. 
Node number of Tra. BC Value of Tra. BC Add c N 
Fixed freedom 
Display input data Generate fle 3DFEM. dat Compute ^S 
4 Element number: © and © 
Display output data Retum Exit NS Node number: 1-12 
Model of three-dimensional program of solid compression 
Display input data: Display output data: 


Fic. 9.38 — Young's modulus and Poisson's ratio in the three-dimensional program of solid 
compression. 
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(4) Step 4 


Figure 9.39 shows one parameter in the fourth line of the parameter input 
window: 

a-coordinate of node: the z-coordinate of each node. The z-coordinate value of 
each node can be input sequentially, and the option ‘Add’ is selected. Then, the 
window ‘Add successfully’ appears, and the option ‘Agree’ is selected to confirm the 
addition. The values entered in this figure were 0 and 0.5. 


a 
Element node is 
Element number in x — || Element number in y |? Element number in = |1 Integration points |8 Material |1 
Young's modulus E. 1.0e6 Poissonsratio v. [0-3 Ko 
x-coordinate ofnode |0 Add bd 
v-coordinate of node Add 
=-coordinate of node Add 
Dis. BC no. P 
Node number of Dis. BC Value of Dis. BC Add o 
Tra. BC no. 
Node number of Tra. BC [=a Add || «c 12N 
Fixed freedom 
| ^M 
Display input data Gener Compute 
| bones BES o | 4 Element number: @ and © 
Display output data Retu Exit XN Node number: 1-12 
PIC:  l— ——— Model of three-dimensional program of solid compression 
Display input data: Display output data: 


Fic. 9.39 — The z-coordinate of node in the three-dimensional program of solid compression. 


(5) Step 5 


Figure 9.40 shows one parameter in the fifth line of the parameter input window: 

y-coordinate of node: The ycoordinate of each node. The ycoordinate value of 
each node can be input sequentially, and the option ‘Add’ is selected. Then, the 
window ‘Add successfully’ appears, and the option ‘Agree’ is selected to confirm the 
addition. The values entered in this figure were 0, 1.5, and 3. 


(6) Step 6 


Figure 9.41 shows one parameter in the sixth line of the parameter input window: 

zcoordinate of node: The zcoordinate of each node. The zcoordinate value of 
each node can be input sequentially, and the option ‘Add’ is selected. Then, the 
window ‘Add successfully’ appears, and the option ‘Agree’ is selected to confirm the 
addition. The values entered in this figure were 0 and —2. 
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By Three-dime 
Element node 

Element number in x 
Young's modulus E 
x-coordinate of node 
y-coordinate of node 
=-coordinate of node 
Dis. BC no. 

Node number of Dis. BC 
Tra. BC no. 

Node number of Tra. BC 
Fixed freedom 


Element number in y | Element number in = |1 Integration points |8 Material |1 
Poisson's ratio v |03 


Element number: © and © 
Node number: 1-12 


Model of three-dimensional program of solid compression 


Fic. 9.40 — The y-coordinate of node in the three-dimensional program of solid compression. 


y-coordinate of node 
2-coordinate of node 
Dis. BC no. 

Node number of Dis. BC 
Tra. BC no. 

Node number of Tra. BC 
Fixed freedom 


[Display input data: 


1 Element number in y |? Element number in = Integration points |8 Material |1 

1.0e6 piak 2810-3 Ag 
0.5 Add M 
3 Add 

o 

L— 


Value of Dis. BC . Adà | 
~~ 
pem 
= 


Element number: © and © 
Node number: 1~12 


Model of three-dimensional program of solid compression 


Display output data: 


Fic. 9.41 — The zcoordinate of node in the three-dimensional program of solid compression. 
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(7) Step 7 


Figure 9.42 shows one parameter in the seventh line of the parameter input 
window: 

Dis. BC no.: the number of displacement boundary conditions. The value entered 
in this figure was 10. 


P3. Three-dimensional program of solid compression - 0 Xx 
Element node 8 
Element number in x — |! Element number in y |? Element number in = |1 Integration points |8 Material |1 
Young s mol 1.0e6 Poksi 10.3 Ag 
x-coordinate ofnode [0-5 Add > 
y-coordinate of node d Add F 10 
z-coordinate ofnode — |? Ada 
Dis. BC no. 10 Sc 
Node number of Dis. BC Value of Dis. BC Add o 
Tra. BC no. 
Node number of Tra. BC Value of Tra. BC Add e DM. 
Fixed freedom 
hc 
Display input data Generate file 3DFEM.dat Compute 
4 Element number: © and © 
Display output data Exit NS Node number: 1-12 
=e Model of three-dimensional program of solid compression 
Display input data: Display output data: 


Fic. 9.42 — Number of displacement boundary conditions in the three-dimensional program 
of solid compression. 


(8) Step 8 


Figure 9.43 shows two parameters in the eighth line of the parameter input 
window: 

Node number of Dis. BC: node number of displacement boundary, specifying the 
node at which a specific displacement occurs. The values entered in this figure were 
1, 2, 3, 4, 7, 8, 9, 10, 11, and 12. 

Value of Dis. BC: values of displacement boundary, specifying the values used for 
displacement boundary in x and y directions, respectively. The value of the dis- 
placement boundary of each specifying node can be input sequentially, and the 
option ‘Add’ is selected. Then, the window ‘Add successfully’ appears, and 
the option ‘Agree’ is selected to confirm the addition. The values entered in this 
figure were 0, 0, and 0. 
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Element node 8 
Element number in x 1 Element number in y |? Element number in = |1 Integration points |8 Material |1 
VELIS 1.0e6 peer OE Ag 
x-coordinate ofnode [0.5 Add = 
y-coordinate ofnode |3 Add 
z-coordiateofnode — |? Add 
Dis. BC no. 10 MX. 
Node number of Dis. BC |1 Value of Dis. BC |0 o o [ au ] o 
Tra. BC no. 
Node number of Tra. BC Vd| Tweed. — O X . Add || ¢ E 
Fixed freedom 
Add successfully M 
Display input data — Compute 
MT NN 4 Element number: © and © 
Display output data l Exit S Node number: 1-12 
E Model of three-dimensional program of solid compression 


Display input data: Display output data: 


Fic. 9.48 — Node number and value of displacement boundary conditions in the 
three-dimensional program of solid compression. 


(9) Step 9 


Figure 9.44 shows one parameter in the ninth line of the parameter input 
window: 

Tra. BC no.: the number of traction boundary conditions. The value entered in 
this figure was 2. 


(10) Step 10 


Figure 9.45 shows two parameters in the tenth line of the parameter input 
window: 

Node number of Tra. BC: node number of traction boundary, specifying the 
node at which specific traction occurs. The values entered in this figure were 5 and 6. 

Value of Tra. BC: values of traction boundary, specifying the values used for 
traction boundary in z-, y-, and z-directions, respectively. The value of the traction 
boundary of each specifying node can be input sequentially, and the option *Add' is 
selected. Then, the window ‘Add successfully’ appears, and the option ‘Agree’ is 
selected to confirm the addition. The value entered in this figure were 0, 0, and 
—1.0e6. 
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E Three-dimensional program of solid compression 


Element node 
Element number in x 


Young's modulus E 
x-coordinate of node 
y-coordinate of node 
z-coordinate of node 


Dis. BC no. 
Node number of Dis. BC o o o Add 
Tra. BC no. 


Node number of Tra. BC Value of Tra BC au | 


Fixed freedom 


Display input data Generate file 3DFEM dat 
Element number: © and © 


Display output data Retum it Node number: 1-12 


Model of three-dimensional program of solid compression 


Display output data: 


Fic. 9.44 — Number of traction boundary condition in the three-dimensional program of solid 
compression. 


f. Three-dimensional program of solid co 
Element node 
Element number in x 


Young's modulus E 
x-coordinate of node 
y-coordinate of node 


z-coordinate of node 
Dis. BC no. 

Node number of Dis. BC 
Tra. BC no. 

Node number of Tra. BC Value of Tra. BC [0 g 
Fixed freedom 


Display input data 


3 Element number: (1) and (2) 
Display output data (i j Node number: 1~12 


Model of three-dimensional program of solid compression 


Display input data: Display output data: 


Fic. 9.45 — Node number and value of traction boundary condition in the three-dimensional 
program of solid compression. 
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Figure 9.46 shows one parameter in the eleventh line of the parameter input 


window: 


Fixed freedom: the number of nodes with the specified displacement. The value 
entered in this figure was 0 by default. 


I3 Three-dimensional program of solid compression 


Element node 8 


Element number in x 1 


"Youngs modulus X 1.0e6 
x-coordinate of node 0.5 Add 
y-coordinate ofnode — |3 ELS 
z-coordinate of node -2 Add 
Dis. BC no. 10 


Node number of Dis. BC |12 Value of Dis. BC |0 o o 
Tra. BC no. 2 
Node number of Tra BC |6 Value of Tra. BC |0 o -1.0e6 
Fixed freedom o 
Display input data Generate file 3DFEM. dat Compute 
Display output data Retum Exit 


4 Element number: © and © 
Bc Node number: 1-12 


Model of three-dimensional program of solid compression 


Display input data: 


Display output data: 


Fic. 9.46 — Fixed freedom in the three-dimensional program of solid compression. 


(12) Step 12 


As shown in figure 9.47, the option ‘Display input data’ is selected, and the 
aforementioned input dataset is summarised on the lower left of the interface. 


(13) Step 13 


As shown in figure 9.48, the option ‘Generate file 3DFEM.dat' is selected, and 
the input dataset is recorded in the file 3DFEM.dat, which is saved in the path of the 
current file. Subsequently, the window 'Generate successfully" appears, and the 
option ‘Agree’ is selected for confirmation. 
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& Three-dimensional program of solid compression - D x 
Element node 8 

Element number in x 1 Element number in y |? Element number in = Integration points |8 Material |1 

Young's modulus E. 1.066 Poisson's ratio» |0-3 Ko 
x-coordinate ofnode [0.5 Add bd 
y-coordinate ofnode — |3 Add 

=-coordinate of node E Add 

Dis. BC no. 10 E 

Node number of Dis. BC. |12 Value of Dis. BC. |0 o o Add o 

Tra. BC no. 2 

Node number of Tra. BC |6 Value of Tra. BC |0 P [1086 Add s 2N 
Fixed freedom jo 


Generate file 3DFEM.dat ‘Compute 
4 Element number: © and © 
Retum Exit NS Node number: 1-12 


Model of three-dimensional program of solid compression 


Display output data: 

12181 
1.0e6 0.3 

5 

153 

2 
10 
1000 2000 3000 4000 
700080009000 10000 
11000 12000 
" 

00-1000000 600 -1000000 


Fic. 9.47 — Input data display in the three-dimensional program of solid compression. 


5 Thise-dlinendoadl progrmobccked — 
Element node 8 
Element number in x 1 Element number in y |? Element number in = |1 Integration points |8 Material |! 
Young's modulus E. 1.0e6 Poisson's ratio v |03 Ko 
eaa of node a 10-5 Add Pd 
y-coordinate ofnode — |3 Ld 
=-coordinate ofnode — |? Add 
Dis. BC no. 10 S 
Node number of Dis. BC |12 Value of Dis. BC |0 o g . Adi | o 
Tra. BC no. 2 
Node number of Tra. BC. |6 Value of Tra. BC |0 o -1.0e6 Add $ 2N 
Fixed freedom o 
Display input data Generate file 3DFEM dat| 
4 Element number: © and © 
Display output data Retum NE Node number: 1-12 
Model of three-dimensional program of solid compression 


100020003000 4000 
7000 80009000 10000 
11000 12000 


5 00 -1000000 600 -1000000 


Fic. 9.48 — Input data generation in the three-dimensional program of solid compression. 
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(14) Step 14 


As shown in figure 9.49, the option ‘Compute’ is selected, and the computa- 
tion program starts to conduct operations using the aforementioned input 
data. Once the computation is completed, the window ‘Compute successfully’ 
appears, and the option ‘Agree’ is selected for confirmation. The output dataset 
above is recorded in the file 3DFEM.res, which was saved in the path of the 
current file. 


8 

Element node 8 

Element number in x 1 Element number in y |? Element number in = |1 Integration points |8 Material |! 

Young s modulus E 1.0e6 Poisson's ratio v [0-3 “Ag 
x-coordinate of node 0.5 Add 

y-coordinate ofnode |3 Aad 

z-coordinate of node r2 Add 

Dis. BC no. 10 RE 

Node number of Dis. BC. |12 Value of Dis. BC |0 o o uL o 

Tra. BC no. 2 

Node number of Tra. BC |6 STeeá. - o x | Hee — Add || ¢ 2N 
Fixed freedom o 

Compute successfully b 
Display input data Compute 
4 Element number: © and © 
Display output data Exit x Node number: 1~12 
Model of three-dimensional program of solid compression 
Display output data: 

12181 

1.0e6 0.3 

0.5 

15 3 
-2 

10 

1000200030004000 

700080009000 10000 

11000 12000 

2 

|5 0 O -1000000 6 0 0 -1000000 


Fic. 9.49 — Computation control in the three-dimensional program of solid compression. 


(15) Step 15 


As shown in figure 9.50, the option ‘Display output data’ is selected, 
and the output data of computed results are displayed on the lower right of the 
interface. 

The option ‘Return’ in the interactive interface, will return to the interface for 
selecting the computation program shown in figure 9.2. The option ‘Exit’ will exit 
the current interactive interface. 
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D Three-dimensional program of solid compression = 0 x 
Element node 8 
Element number in x 1 Element number in y |? Element number in = |1 Integration points |8 Material |! 
Young's modulus E 1.0e6 Poisson's ratio y |0-3 Koa 
x-coordinate ofnode — |05 Add b d 
v-coordmateofnode |3 Add 
z-coordinate of node r2 Add 
Dis. BC no. 10 S 
Node number of Dis. BC. |12 Value of Dis. BC. |0 o o [oe | o 
Tra. BC no. 2 
Node number of Tra. BC |6 Value of Tra. BC |0 0 [ross (pad o 2N 
Fixed freedom o 


Element number: © and © 


4 

Í Display output data | Retum Exit bs Node number: 1-12 

Model of three-dimensional program of solid compression 
p ^ 
R2161 Node xdisp y-disp  z-disp 
1.0e6 0.3 1 0.0000E+00 0.0000E+00 0.0000E+00 
jo .5 2 0.0000E+00 0.0000E+00 0.0000E«00 
p 15 3 3 0.0000E+00 0.0000E+00 0.0000E+00 
p -2 4 0.0000E«00 0.0000E+00 0.0000E+00 
10 5 -0.3438E+00 0.3035E-16 -0.4332E+01 
1000 2000 3000 4000 6 0.3438E+00 0.1720E-16 -0.4332E+01 
7000 80009000 10000 7 0.0000E+00 0.0000E+00 0.0000E+00 
11000 12000 8 0.0000E+00 0.0000E+00 0.0000E+00 
" 9 0,0000E+00 0.0000E+00 0.0000E+00 
5 00 -1000000 6 00 -1000000 10 0.0000E+00 0.0000E+00 0.0000E+00 y 


Fic. 9.50 — Output data display in the three-dimensional program of solid compression. 


9.4 Exercises 


9.4.1 Using the one-, two-, and three-dimensional programs provided above, com- 
pute the displacements and stresses for the following one-, two-, and 
three-dimensional models with various values of the parameters a, b, and c, as shown 
in figure 9.51. 

9.4.2 Using the analytical results of the theoretical models, verify the correctness of 
the programs in this chapter, and discuss the influence of different mesh densities on 
the accuracy of the solutions. 

9.4.3 Based on the one-dimensional program of beam deformation, use shape 
functions, and numerical integration to compute the results. 

9.4.4 Based on the two-dimensional program of plane strain problem, compute the 
solutions for the plane stress model. 
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(a) One-dimensional model 


F=-1x10° | 


(c) Three-dimensional model 


Fic. 9.51 — One-, two-, and three-dimensional models with variable geometric sizes (a, b, 
and c). 


Appendix A. Keyword Index 


The following is a list of frequently used keywords in this book, which can be quickly 
indexed to the specified chapters. It should be noted that some keywords appear and 
are used in multiple chapters, which will be given simultaneously. 


Chapter 1 


Computational mechanics 
Numerical algorithm 
Continuous system 

Finite element method 


Element 


Chapter 2 


Three-dimensional problem 
Plane stress problem 
Axisymmetric problem 
Strain 

Geometric equations 


Equilibrium equations 


Chapter 3 


Strong form 
Differential 


Arbitrary function 
Derivative 

Shape function 
Load 


Numerical method 
Numerical model 
Discrete system 

Finite element analysis 
Mesh 


Two-dimensional problem 
Plane strain problem 
Displacement 

Stress 

Constitutive equations 


Boundary conditions 


Weak form 
Integration 


Integration by parts 
Galerkin method 
Stiffness 


Element location vector 
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Chapter 4 


Lagrange element 
Triangle 

Triangle element 
Rectangle 
Rectangle element 
Tetrahedron 


Tetrahedron element 


Chapter 5 


Isoparametric mapping 
Isoparametric element 
Global coordinate 
Local coordinate 


Derivative 


Chapter 6 


Elasticity problem 
Two-dimensional problem 
Strain 

Geometric equations 
Equilibrium equations 


Weak form 


Chapter 7 


Algebraic equations 
Sti 
Decomposition 
Row 


Column 


ness matrix 


E 


Hexagon 

Hexagon element 
Area coordinate 
Volume coordinate 
Linear element 
Quadratic element 


Cubic element 


Jacobian matrix 
Numerical integration 
Gaussian quadrature 
Integration point 
Weight 


'Three-dimensional problem 
Displacement 

Stress 

Constitutive equations 


Boundary conditions 


Stiffness equations 


Upper triangular matrix 
Lower triangular matrix 
Gaussian elimination 
Forward elimination 


Back substitution 
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Chapter 8 


Error estimation 
Adaptive finite element method 


Error tolerance 
Chapter 9 


Program 
Geometric model 


Arrays 


Superconvergent solutions 
Superconvergent patch recovery 


Solution refinement 


Modules 
Input data 


Output solutions 
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Appendix B. Matrix Calculation 


A feature of the finite element method is that systems of equations are compactly 
expressed in matrix form, and thus matrix operations can be used to simplify the 
calculations. The basic matrix concepts and operations used in this book include 
matrix definitions, addition or subtraction, transpose of a matrix, transpose of a 
product, inverse, symmetric matrices, and partitioning of matrices. 


B.1 Definition 


The following linear relations between a set of variables x; and b; (i = 1, 2,3; j = 1, 2) 


atı + a212 + a4323 = by 


(B.1) 


213 + 22% + 02333 = b2 


can be written as: 


Alz} ={b} (B.2) 


or 


Ax=b (B.3) 


asjaga e 3] po yada bbe (y= : (B.A) 
| | Us) 


5j, 422 423 d 


'The above notation involves the concept of a matrix and the operation of matrix 
multiplication. Matrices are defined as ‘arrays of number’ of the type in equation 
(B.4). Matrices containing a single column of numbers are often referred to as vector 
or column matrices, whereas a matrix with multiple columns and rows is called a 
rectangular matrix. The multiplication of a matrix by a column vector, as in 
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equation (B.2), is defined by the equivalence of the left and right sides of equations 
(B.1) and (B.2). Boldface characters are used to denote both vectors and matrices. 

If another relationship to equation (B.1), using the same a, but a different set of 
xz and b, we have the following linear relations: 


aay == 0125, + 01325 = b 
/ / / l (B.5) 
0517, F 2213 + 02334 = b 
then equations (B.1) and (B.5) can be simultaneously expressed as: 
[AX] = [B] or AX—B Ba 
where 
Ti Ti 
a a a 
A=[4] zi Ks a| x= =| ay ot |: 
ids T; 0x (B.7) 
b, b 
B = [B] -| : 1 
b b 
or equivalently 
autt. at + 2 NER E E bi b 
be TE ay 0, + i -=B =[8]=| b, (B.8) 


Two matrices are equal if and only if their corresponding entries are equal. 

The multiplication of full matrices is defined in equation (B.6), and evidently, it is 
meaningful only if the number of columns in A is equal to the number of rows in 
X. One property that distinguishes matrix multiplication is that, in general, 


AX + XA 


That is, unlike scalar multiplication, matrix multiplication is not commutative. 


B.2 Matrix Addition or Subtraction 


If relations of the form from equations (B.1) and (B.5) are added, then we have: 


H a13 (a3 4 ) =b + b 


alt T xi) + di»(25 4 Fas 
ie 93 (33 ar zh) = b» a b, 


ayy (T1 T i) sir a2 (m 1 


(B.9) 


+ 25) 4 
+ 25) 4 


which can also be written as: 
Ax -- Ax! =b +b' 


if we define the addition of two matrices as the matrix, each entry of which is the 
sum of the corresponding entries of the two matrices. Clearly, this can be done only if 
the matrices have the same dimension (same number of rows and columns). For 
example, 
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Qj] a2 bi b a + by a2 + bye ĉii C12 
anı az | + | ba, b2 | = | az + b2 a22 + 023 | = | €i. €22 
a31 032 bs, b32 d31 + b31 a32 + b32 C31 C32 
Or 
A+B=C (B.10) 


implies that every entry of C is equal to the sum of the corresponding entries of 
A and B. Matrix subtraction is defined similarly. 


B.3 Transpose 


'The transpose of a matrix is obtained by interchanging the rows and columns of the 
given matrix and is denoted by the superscript T. For example, 


T a1 agi 
Qjj a2 3 
= 012 a22 (B.11) 
a21 022 423 
013 23 


In mechanics problems, we often encounter a number of quantities, such as force, 
which can be expressed in vector form: 


h 
f= ; (B.12) 
hh 
This, in turn, is often associated with the same number of displacements given by 
another vector: 


u=} . (B.13) 
Un 


It is known that work is represented as the sum of the products of force and 
displacement: 


S-Y hu 
k=1 


Clearly, the transpose is useful here, as, by the rule of matrix multiplication, we 
have 


UL 
S= {A b hy" s —flu-ulf (B.14) 


Un 
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B.4 Transpose of a Product 


An operation that sometimes occurs is taking the transpose of a matrix product. 
The following equation can be proved by examining the corresponding entries of the 
matrices on the two sides: 


(AB)! = BT AT (B.15) 
B.5 Inverse 


If in equation (B.3), the matrix A is ‘square’, that is, it represents the coefficients of 
simultaneous equation (B.1), the number of which is equal to the number of 
unknowns z, then, under certain assumptions, it is possible to solve for unknowns in 
terms of the known coefficients b. This solution can be written as: 


x—-A b (B.16) 


where matrix A‘ is known as the ‘inverse’ of the square matrix B. Clearly, A‘ is 
also square and of the same dimension as B. 
We can obtain equation (B.12) by multiplying both sides of equation (B.3) by 
A* 
A A-I-AA! (B.17) 
where Iis the ‘identity’ matrix, that is, it has 0 at all off-diagonal positions and 1 on 
the diagonal. 


B.6 Symmetric Matrices 


Symmetric matrices are often encountered in structural problems. If the general 
entry of a matrix A is denoted by aj, then A is called symmetric if 


aj =a; or A-—AT. (B.18) 
A symmetric matrix must be square. It can be shown that the inverse of a 


symmetric matrix is also symmetric: 


A™= (A) = AT. (B.19) 


B.7 Partitioning 


The product of two matrices A and B can also be expressed in terms of their blocks. 
For instance, if 
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bu bie 
a1 a2 013 0414 Q5 b b22 
A= |a 09 a3 ax azg |, B= | bs, b» (B.20) 
a31 a32 | (33 a34 a35 ba ba 
b51 b52 


then A and B can be multiplied by partitioning each matrix into submatrices, 
indicated by the lines, and applying the rules of matrix multiplication first to each of 
these submatrices, as if it were a scalar, and then carrying out further multiplication 
in the usual manner. Specifically, A and B can be written as: 


An Arm Bi 
= B= B.21 
ie Az B, pod 
and then, the product A B can be expressed as: 


(B.22) 


AB = co Mur 


A»Bi T Az3B» 


'This product can be verified as representing the complete product by further 
multiplication. 

The essential feature of partitioning is that the size of subdivisions has to be such 
as to make the products of the type A,,;B, meaningful, i.e., the number of columns in 
Aj, must be equal to the number of rows in B,, etc. If the above definition holds, 
then all further operations can be conducted on partitioned matrices, treating each 
partition as if it were a scalar. 

It should be noted that any matrix can be multiplied by a scalar (number). Here, 
obviously, the requirements of equality of appropriate rows and columns no longer 
apply: 


Aa Aay Aas Aaya Aas 
AA = 4A = Aan, A22 À.d23 À.024 A025 (B.23) 
azı azz Ad33 Adg4  Ad35 


If a symmetric matrix is divided into an equal number of submatrices Aj; in rows 
and columns, then we have 


Aj =A; (B.24) 


Appendix C. Summary of Elements 
and Shape Functions 


The following are the elements and shape functions for the one-, two-, and 
three-dimensional cases introduced in this book. For the convenience of the reader, 
they are summarised to be indexed and compared. 


C.1 One-Dimensional Lagrange Element 


Linear element with two nodes 
J J 


N(x) 21 = and M(x) = - with x wx and he= ay — zf (C1) 


€ € 


Quadratic Lagrange element 


MO) = 80) = Ca 23-0 
MO = IO = REET Sete) (c3) 
MG = BO = EER = 1-8 


C.2 Two-Dimensional Triangle Element 


Triangle with three nodes 


1 
N(x, y) = zaet baz + cay), a=1,2,3 
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% = RY- BY, b=y-y, =B- H 
with a? = By — 31, b2 =Y — Y, & = Tt — 3 (C.3) 
a3 = UY — MY, bg = YY, B=- Ü f 
and A= (aby + 2909 + 23 03)/2 


Quadratic triangle element 
Corner nodes: 


N, = (2La — 1)La, a=1, 2,3 (C.4) 
Mid-side nodes: 
Ny = 4LilLo, N; =4IoL3, Neg = 403l, (C.5) 
Cubic triangle element 
Corner nodes: 
Ny = TOL —1)(3L,—2)L, a=1, 2,3 (C.6) 


Mid-side nodes: 
9 9 9 
M=5hb(8h-1), NM -5lL(L-1), No=5lal(3l2-1) 
(C.7) 
9 9 9 
N; = z I2la(3Ls —1), M= z BBL —1), M= z BABL - 1) 


Internal node: 


Nio = 27 L Lela (C.8) 


C.3 Two-Dimensional Rectangle Element 


Linear rectangle element with four nodes 


1 1 
Nı —M EU M 2410-90 - n) 
1 1 (C.9) 
M-i0*O0e5, M-10-000 
Quadratic rectangle element 
Corner nodes: 
N,- iE (E+E) MHN), @=1, 2, 3,4 (C.10) 


(4 


Appendix C. Summary of Elements and Shape Functions 177 


Mid-side nodes: 


1 
Ca = 9, Ny 3n (1 & (n 4 1); a=5,7 
1 (C.11) 
n, = 0, Na= 5 (E+€,) (1-5), a=6, 8 
Center node: 
Ny = (1 - 8) (1- P) (C.12) 
C.4 Three-Dimensional Tetrahedron Element 
Linear tetrahedron element with four nodes 
1 
N.(2, yY, z) = ey T bat + cay + daz), a=1, 2, 3, 4 (C.13) 
Quadratic tetrahedron element 
Corner nodes: 
N, = La (2La— 1), a=1, 2, 3, 4 (C.14) 
Mid-side nodes: 
N; = 42l, Ne = 4L;Li, N= 4L4L (C.15) 
Ng = 4LəL3, No = 4LzL4, Nin = 4L4L2 ` 
Cubic tetrahedron element 
Corner nodes: 
1 
Nai = 5 (3L, — 1) (3La— 2), a=1, 2, 3,4 (C.16) 
Mid-side nodes: 
9 9 9 
N; = 5I bL —1, N= z Ilo -1) N= z I las —1) 
9 9 9 
Ng = z bL -1, N= z Ils -1) Mo= z Peds (GL - 1) 
(C.17) 


9 9 9 
Ny = z Ils (3La =1); Ne= z Ila (3L. -1, M= z 234a (3L - 1) 


9 9 9 
Na = z Lla (3L —1, M5; = z Isla Ls —1) Ne= z I lala — 1) 
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Mid-face nodes: 
Niz = 27 Li IoL3, Nig = 27 L4 Li Lo 


M= Misli. Noy — 2Tolabs oe 
C.5 Three-Dimensional Hexahedron Element 
Hexahedron with eight nodes 
M-$0-80ü-90-0. M-;0-800-9 49) 
M$-$iü-)ü-59ü-0, M-;0«00-20«0 
A : (C.19) 
Asc. 5 Ita) =O, Nr-sQ-5 Uem (EO 
M= -AUU M= A-E (+r) 040 
Quadratic hexahedron element 
Corner nodes: 
Ny = nt (E+E) 0) (8G). G51, 2-4 8 (C.20) 


Mid-side nodes: 


Ca = 0, N, = nt (1—&) (nm) (C4+6.), a- 9, 10, 11, 12 
no =0, Nasi (E) (Lo) CHi), a= 13, 14, 15,16 (C21) 
Ca = 0, No Zen (E+E) (n+) (I=L), a 17, 18, 19, 20 
Mid-face nodes: 
Ča = Na =Q, N= a-&8)0-45) (Cti), a=21, 22 
Na = a = 0, Na = (E+ éa) (1-7?) 1-0), a= 23, 24 (C.22) 
Ca = ča = 0, Na 5 (1- 8) (nn) 1-@), a= 25, 26 


Center node: 


Nv=(1-&@) (l-A (1-0), a=27 (C.23) 


Appendix D. Gaussian Integration 
Points and Weights 


For one-dimensional problems, Gaussian integration is given by: 


[rows Ere (D.1) 


with an appropriate choice of the integration points č; and weights wj. The Gaussian 
quadrature formula in equation (D.1) using n integration points is exact in the case 
of polynomials of degree at most 2n — 1. 

For high-dimensional problems, the number of integration points in each 
dimension abides by the above rules. Accordingly, for two-dimensional problems, we 
can integrate into two directions using 


a [ Pii (é, n) Jel (é, " )dédn ~ ow Y Y (oh) Je [usi Nn) Um Ur, (D.2) 


n=l m=1 


where m denotes quadrature points in the c direction and n quadrature points in the 
y direction. For three-dimensional problems, we can integrate into three directions 
using 


4 l / l ri Ga Qin Sagara 


LN M (D.3) 
22399 f( Ems Nn CII (Ems Nn, $1) Wm Wn 


l=1 n=l m=1 


where m denotes quadrature points in the c direction, n quadrature points in the 7 
direction, and l quadrature points in the ¢ direction. 

For complex integrands expressed by higher-order polynomials, more integration 
points are required to obtain exact formulas. The following decimal forms of these 
integration points and weights are provided. 
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n=1 


Integration points 
0.0000000000000000 


n=2 


Integration points 
~0.5773502691896250 
0.5773502691896250 


n=3 


Integration points 
—0.7745966692414830 
0.0000000000000000 
0.7745966692414830 


n=4 


Integration points 
-0.8611363115940520 
-0.3399810435848560 
0.3399810435848560 
0.8611363115940520 


n=5 


Integration points 
—0.9061798459386640 
-0.5384693101056830 


0.0000000000000000 
0.5384693101056830 


0.9061798459386640 


Weights 


2.0000000000000000 


Weights 
1.0000000000000000 
1.0000000000000000 


Weights 
0.5555555555555550 
0.8888888888888880 
0.5555555555555550 


Weights 
0.3478548451374530 
0.6521451548625460 
0.6521451548625460 
0.3478548451374530 


Weights 
0.2369268850561890 
0.4786286704993660 


0.5688888888888880 
0.4786286704993660 


0.2369268850561890 


n=6 


Integration points 
-0.9324695 142031520 
-0.6612093864662640 
-0.2386191860831960 
0.2386191860831960 
0.6612093864662640 
0.9324695 142031520 


n=7 


Integration points 
~0.9491079123427580 
-0.7415311855993940 
-0.4058451513773970 
0.0000000000000000 
0.4058451513773970 
0.7415311855993940 
0.9491079123427580 


n=8 


Integration points 
-0.9602898564975360 
-0.7966664774136260 

—0.5255324099163290 
-0.1834346424956490 
0.1834346424956490 
0.5255324099163290 
0.7966664774136260 
0.9602898564975360 
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Weights 
0.1713244923791700 
0.3607615730481380 
0.4679139345726910 
0.4679139345726910 
0.3607615730481380 
0.1713244923791700 


Weights 
0.1294849661688690 
0.2797053914892760 
0.3818300505051180 
0.4179591836734690 
0.3818300505051180 
0.2797053914892760 
0.1294849661688690 


Weights 
0.1012285362903760 
0.2223810344533740 
0.3137066458778870 
0.3626837833783620 
0.3626837833783620 
0.3137066458778870 
0.2223810344533740 
0.1012285362903760 
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n=9 


Integration points 
~0.968 1602395076260 
-0.8360311073266350 
-0.6133714327005900 
-0.3242534234038080 
0.0000000000000000 
0.3242534234038080 
0.6133714327005900 
0.8360311073266350 
0.9681602395076260 


Weights 
0.0812743883615744 
0.1806481606948570 
0.2606106964029350 
0.3123470770400020 
0.3302393550012590 
0.3123470770400020 
0.2606106964029350 
0.1806481606948570 
0.0812743883615744 


n= 10 


Integration points 
-0.9739065285171710 
-0.8650633666889840 
—0.6794095682990240 
-0.4333953941292470 
-0.1488743389816310 
0.1488743389816310 
0.4333953941292470 
0.6794095682990240 
0.8650633666889840 
0.9739065285171710 


Weights 
0.0666713443086881 
0.1494513491505800 
0.2190863625159820 
0.2692667193099960 
0.2955242247147520 
0.2955242247147520 
0.2692667193099960 
0.2190863625159820 
0.1494513491505800 
0.066671344308688 1 


n= 11 


Integration points 
-0.9782286581460570 
—0.8870625997680950 
-0.7301520055740490 
-0.5190961292068110 
-0.2695431559523440 
0.0000000000000000 
0.2695431559523440 
0.5190961292068110 
0.7301520055740490 
0.8870625997680950 
0.9782286581460570 


Weights 
0.0556685671161736 
0.1255803694649040 
0.1862902109277340 
0.2331937645919900 
0.2628045445102460 
0.2729250867779000 
0.2628045445102460 
0.2331937645919900 
0.1862902109277340 
0.1255803694649040 
0.0556685671161736 
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n= 12 


Integration points 
-0.9815606342467190 
-0.9041172563704740 
-0.7699026741943040 
-0.5873179542866170 
-0.3678314989981800 
-0.1252334085114680 
0.1252334085114680 
0.3678314989981800 
0.5873179542866170 
0.7699026741943040 
0.9041172563704740 
0.9815606342467190 


Weights 
0.0471753363865118 
0.1069393259953180 
0.1600783285433460 
0.2031674267230650 
0.2334925365383540 
0.2491470458134020 
0.2491470458134020 
0.2334925365383540 
0.2031674267230650 
0.1600783285433460 
0.1069393259953180 
0.0471753363865118 


n= 13 


Integration points 
-0.9841830547185880 
-0.9175983992229 710 
-0.8015780907333090 
-0.6423493394403400 
-0.4484927510364460 
—0.2304583159551340 
0.0000000000000000 
0.2304583159551340 
0.4484927510364460 
0.6423493394403400 
0.8015780907333090 
0.9175983992229770 
0.9841830547185880 


Weights 
0.0404840047653158 
0.0921214998377284 
0.1388735102197870 
0.1781459807619450 
0.2078160475368880 
0.2262831802628970 
0.2325515532308730 
0.2262831802628970 
0.2078160475368880 
0.1781459807619450 
0.1388735102197870 
0.0921214998377284 
0.0404840047653158 
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n= 14 


Integration points 
-0.9862838086968120 
-0.9284348836635 730 
~0.8272013150697650 
-0.6872929048116850 
—0.5152486363581540 
~0.3191123689278890 
-0.1080549487073430 
0.1080549487073430 
0.3191123689278890 
0.5152486363581540 
0.6872929048116850 
0.8272013150697650 
0.9284348836635730 
0.9862838086968120 


Weights 
0.0351194603317518 
0.0801580871597602 
0.1215185706879030 
0.1572031671581930 
0.1855383974779370 
0.2051984637212950 
0.2152638534631570 
0.2152638534631570 
0.2051984637212950 
0.1855383974779370 
0.1572031671581930 
0.1215185706879030 
0.0801580871597602 
0.0351194603317518 


n=15 


Integration points 
—0.9879925180204850 
-0.9372733924007060 
-0.8482065834104270 
~0.7244177313601700 
~0.5709721726085380 
-0.3941513470775630 
-0.2011940939974340 
0.0000000000000000 
0.2011940939974340 
0.3941513470775630 
0.5709721726085380 
0.7244177313601700 


0.8482065834104270 
0.9372733924007060 


0.9879925 180204850 


Weights 
0.0307532419961172 
0.0703660474881081 
0.1071592204671710 
0.1395706779261540 
0.1662692058169930 
0.1861610000155620 
0.1984314853271110 
0.2025782419255610 
0.1984314853271110 
0.1861610000155620 
0.1662692058169930 
0.1395706779261540 


0.1071592204671710 
0.0703660474881081 


0.0307532419961172 


Appendix D. Gaussian Integration Points and Weights 


n= 16 


Integration points 
-0.9894009349916490 
~0.9445750230732320 
—0.8656312023878310 

-0.7554044083550030 

-0.6178762444026430 
~0.4580167776572270 
~0.2816035507792580 

-0.0950125098376374 
0.0950125098376374 
0.2816035507792580 
0.4580167776572270 
0.6178762444026430 
0.7554044083550030 
0.8656312023878310 
0.9445750230732320 
0.9894009349916490 


Weights 
0.0271524594117540 
0.0622535239386478 
0.0951585116824927 
0.1246289712555330 
0.1495959888165760 
0.1691565193950020 
0.1826034150449230 
0.1894506 104550680 
0.1894506104550680 
0.1826034150449230 
0.1691565193950020 
0.1495959888165760 
0.1246289712555330 
0.0951585116824927 
0.0622535239386478 
0.0271524594117540 


n= 17 


Integration points 
-0.9905754753144170 
-0.9506755217687670 

-0.8802391537269850 
-0.7815140038968010 
-0.6576711592166900 
-0.5126905370864760 
-0.3512317634538760 
-0.1784841814958470 
0.0000000000000000 
0.1784841814958470 
0.3512317634538760 


0.5126905370864760 
0.6576711592166900 


0.7815140038968010 
0.8802391537269850 
0.9506755217687670 
0.9905754753144170 


Weights 
0.0241483028685479 
0.0554595293739872 
0.0850361483171791 
0.1118838471934030 
0.1351363684685250 
0.1540457610768100 
0.1680041021564500 
0.1765627053669920 
0.1794464703562060 
0.1765627053669920 
0.1680041021564500 


0.1540457610768100 
0.1351363684685250 


0.1118838471934030 
0.0850361483171791 
0.0554595293739872 
0.0241483028685479 
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n=18 


Integration points 
-0.9915651684209300 
-0.95582394957139770 

0.8926024664975550 
-0.8037049589725230 
-0.6916870430603530 
—0.559777083107394'70 
-0.4117511614628420 
-0.2518862256915050 
~0.0847750130417353 
0.0847750130417353 
0.2518862256915050 
0.4117511614628420 
0.55977083 10739470 
0.6916870430603530 
0.8037049589725230 
0.8926024664975550 
0.9558239495713970 
0.9915651684209300 


Weights 
0.0216160135264833 
0.0497145488949698 
0.0764257302548890 
0.1009420441062870 
0.1225552067114780 
0.1406429146706500 
0.1546846751262650 
0.1642764837458320 
0.1691423829631430 
0.1691423829631430 
0.1642764837458320 
0.1546846751262650 
0.1406429146706500 
0.1225552067114780 
0.1009420441062870 
0.0764257302548890 
0.0497145488949698 
0.0216160135264833 


n= 19 


Integration points 
-0.9924068438435840 
-0.9602081521348300 

-0.9031559036148170 
—0.8227146565371420 
-0.72096617773352290 

-0.6005453046616810 
-0.4645707413759600 

-0.3165640999636290 
—0.1603586456402250 
0.0000000000000000 
0.1603586456402250 
0.3165640999636290 
0.4645707413759600 
0.6005453046616810 
0.7209661773352290 
0.8227146565371420 
0.9031559036148170 
0.9602081521348300 
0.9924068438435840 


Weights 
0.0194617882297264 
0.0448142267656996 
0.0690445427376412 
0.0914900216224500 
0.1115666455473330 
0.1287539625393360 
0.1426067021736060 
0.1527660420658590 
0.1589688433939540 
0.1610544498487830 
0.1589688433939540 
0.1527660420658590 
0.1426067021736060 
0.1287539625393360 
0.1115666455473330 
0.0914900216224500 
0.0690445427376412 
0.0448142267656996 
0.0194617882297264 
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n = 20 


Integration points 
-0.9931285991850940 
-0.9639719272779130 
-0.9122344282513260 
-0.8391169718222180 
-0.7463319064601500 
-0.6360536807265150 
~0.5108670019508270 
~0.3737060887154190 
-0.22777858511416450 
-0.0765265211334973 
0.0765265211334973 
0.2277858511416450 
0.3737060887154190 
0.5108670019508270 
0.6360536807265150 
0.7463319064601500 
0.8391169718222180 
0.9122344282513260 
0.9639719272779130 
0.9931285991850940 


Weights 
0.0176140071391521 
0.0406014298003869 
0.0626720483341090 
0.0832767415767047 
0.1019301198172400 
0.1181945319615180 
0.1316886384491760 
0.1420961093183820 
0.1491729864726030 
0.1527533871307250 
0.1527533871307250 
0.1491729864726030 
0.1420961093183820 
0.1316886384491760 
0.1181945319615180 
0.1019301198172400 
0.0832767415767047 
0.0626720483341090 
0.0406014298003869 
0.0176140071391521 
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Chapter 1 


1.4.1 Describe the treatment of ‘standard discrete problems’. 

Solutions: 

The existence of a unified treatment of ‘standard discrete problems’ leads us to the 
first definition of the finite element process as a method of approximation to con- 
tinuum problems such that: 


(a) 
(b) 


the continuum is divided into a finite number of parts (elements), the behaviour 
of which is specified by a finite number of parameters, and 

the solution for the complete system as an assembly of its elements precisely 
follows the same rules as those applicable to standard discrete problems. 


1.4.2 Summarise the computation procedure of the finite element method. 
Solutions: 

It should be clear to the reader that the finite element solution of a problem always 
follows a standard methodology. The solution process for any problem type is always 
performed by the following steps: 


(1) 


Define the problem to be solved in terms of differential equations. Construct the 
integral form for the problem as virtual work, variational, or weak formulation. 
Select the type and order of the finite elements to be used in the analysis. 
Define the mesh for the problem. This involves the description of the node and 
element layout, as well as the specification of boundary conditions and 
parameters for the formulation used. 

Compute and assemble the element arrays. The particular virtual work, vari- 
ational, or weak form provides the basis for computing the specific relationships 
of each element. 

Solve the resulting set of linear algebraic equations for the unknown parame- 
ters. For static problems, this requires the solution of linear algebraic equations. 
Output the results for the nodal and element variables. Graphical outputs are 
also useful in this step. 
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Chapter 2 


2.8.1 Summarise the basic variables and equations of three- and two-dimensional 
elasticity. 

Solutions: 

1. Three-dimensional elasticity 

Basic variables: 


(1) Displacement components: u, v, w 


(2) Strain components: Er, £y, €z, Vays Vyzr Vow 


(3) Stress components: Os, Oy, Oz, Tay, Tyz, Tex 
Basic equations: 


(1) Geometric equations: 


o 
T 
^ 0 5, 9 
êy 0 0 S fu 
& l Oz 
e= " = ð fa) v 
Yzy cu. Ao ak 0 w 
Viz Oy Or 
Oz Oy 
o O 
ir | cem 
Oz Ox 
(2) Constitutive equations: 
o = De 
(1— v) v U 0 0 0 
v (1 — v) U 0 0 0 
: E v v (1— v) 0 0 0 
imd MES 0 0  (1-20/2 0 0 
0 0 0 0 (1— 2v)/2 0 
0 0 0 0 0 (1— 2v)/2 
(3) Equilibrium equations: 
O0; | Oty; Ota 
boss] 
Ox Oy 2 Oz P 
Ota Ooy Oty u 
Ox Oy F Oz a 
OTr: Oty, OG, is cep 
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2. Two-dimensional elasticity (Plane stress and stain case) 
Basic variables: 


(1) Displacement components: u, v 
(2) Strain components: &;, &y, êz, y, 
(3) Stress components: Or, Oy, Oz, tay 


Basic equations: 


(1) Geometric equations: 


ð 
0 
Ox 
Er 0 0 
gm | = oy HL : = Syu T £; 
T p ete Ts 
Pay 0 
a 
Oy Or 
(2) Constitutive equations: 
Or 1 v 0 0 Er 
Jol. E v 1 0 0 by 
Plane stress: o | (i-3j|0 0 0 0 5 
ay 0 0 0 (1—9)2] lta 
Or (1— v) v v 0 Ez 
_ Joy V E v (1 — v) v 0 £y 
Plane strain: d a i » (1— v) 0 hs 
Tay 0 0 0 (1— 20)/2 | UY 
(3) Equilibrium equations: 
Oo, | Oye 
+b, =0 
Ox i Oy 
Ot, Oo, 
—— L——rb,-0 
Ox T Oy " 


2.8.2 Using the formulation presented in this chapter, provide the basic variables 
and equations of axial deformation for the one-dimensional elastic cantilever beam 
shown in figure 2.5. The coordinates of the two ends are z — a and z — b, and the 
axial uniform load is b, 

Solutions: 

1. Basic variables: 


(1) Displacement component: ur 
(2) Strain component: £, 
(3) Stress component: 6; 
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Basic equations: 


; , Qu 
(1) Geometric equation: ¢, = — 


Ox 


(2) Constitutive equation: e; = Ee; 


(3) Equilibrium differential equation: oe +b: =0 
z 


Chapter 3 


3.6.1 Describe the implementation steps for the weak form of equivalent integration 
Solutions: 

A weak form for a set of differential equations is obtained through the following 
steps: 


(1) Multiply each equation by an appropriate arbitrary function. 

(2) Integrate this product over the space domain of the problem. 

(3) Use integration by parts to reduce the order of the derivatives to a minimum. 
(4) Introduce boundary conditions if possible. 


3.6.2 Using the Galerkin method, for example 3.1 with E(x) = E - x, bz(x) = bz: x, 
compute the global stiffness matrix K and load matrix F. 

Solutions: 

The conditions of these problems are 


10x for0cz«5 
0 for5b«r«10 
Boundary conditions: u(0) 2 0. and £,(10) = 25 
10 10 
[;) f2] 
T hr f w(z)b,zdx — w(10)25 = 0 


a EX 
o Ox Or 0 


Loading: b; = { 


Weak form: 


Consider an approximate solution: 


The weak form is 


N 10 10 
E pNn-ly,gqNm-1 yNm 
nade = (=) bxd + 25 
x um) (Gu € i T EE 


The coefficients of stiffness matrix and load matrix are 


10 
E qzNmTn-2 Emn 
"nem — TAR (=) dz = 
n 100 "1g T mcn 
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5 gum 1000 /1\”+? 
n >) 1 25 = = 2 
f n (=) Orda + 25 s) 4-25 


For example, m — n — 2, E — 1000, b, — 10 the stiffness matrix and load matrix 
are given by: 


200 
K= ku k| _ | 500 667 j= hl Ja 
~ [ko ka| [667 10007 \f ) 325 
8 


3.6.3 Using the finite element method, for example 3.2 with E(z) = E-z, 


ba(x) = bz - x, compute the global stiffness matrix K and load matrix f. 
Solutions: 


We divided the domain into four equal finite elements, so the length of the element is 


he = £5 — tf = 2.5, =r- r? 191,9,3,4 


: a! 
The shape functions are M(x’) = 1 — = M(x’) = AD 
The derivatives of shape functions are 
dN, dN, 1 1 0.4, dN. dN, 1 1 0.4 
dr da’ he 2.5 dr dz he 25 
'The element stiffness matrix and load matrix are 
dN, 
Ne TT 2.5 —0.4 
[ke = i da’ Bs] oN ar = n { Jes-04 0.4 da 
o |2N da! da 0 0.4 
dz' 


zs 0.16 0.16 ; 
= f Er dx 
0 0.16 0.16 


(1) Element stiffness matrix 


Element 1: 


2.5 
1 [016 016],, E[1 -1] [500 -500 
i -f en ee 08 | =3[4, 1 | F e 500 | 
Element 2: 


i - [^ nwalan oue - | : Fl 
0 


0.16 0.16 


1500  —1500 
—1500 1500 


194 Appendix E. Exercise Solutions 


Element 3: 
2.5 
3 j 0.16 0.16 , 9E 1 -1| | 2500 —2500 
K =f Ec 5) [oae 0.16 ism 2 |—1 1| |-2500 2500 
Element 4: 
2.5 
4 7 0.16 0.16 ees 1 -1| | 3500 -—3500 
K =f EG! 1) | 916 on ee 2 |—1 1 | |-3500 3500 
(2) Load matrix 
Element 1: 
25| 1 d 
= f 25 hast 10.42 
f =| a haidi s 
2.5 
Element 2: 
25 | 1 d 
2 25 F , [41.67 
f -f a b, (a! + 2.5) dx SA 
.5 
Element 3: 
J 
5 2.5 let 0 
a / EN 
f = FACED s 
2.5 
Element 4: 


(3) Global stiffness matrix and load matrix 
Element location vector: 
at ={1 2}, ={2 3}, 4 =1{3 4}, à ={4 57 


The element stiffness matrix and load matrix can be integrated into the global 
stiffness matrix and load matrix: 
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1 2 3 4 5 
1| 500 —500 0 0 0 
2 | —500 500+ 1500 —1500 0 0 
K=3 0 —1500 1500 + 2500 —2500 0 
4 0 0 —2500 25004-3500  —3500 
5 0 0 0 —3500 3500 
500 | —500 0 0 0 
—500 2000  —1500 0 0 
= 0 —1500 4000  —2500 0 
0 0 —2500 6000  —3500 
0 0 0 —3500 3500 
1 10.42 10.42 
2 | 20.83 + 41.67 62.5 
f=3 52.08 + 0 = 4 52.08 
A 0+0 0 
5 0 0 


3.6.4 Summarise the functions of the element location vector 
Solutions: Omission. 


Chapter 4 


4.6.1 One-dimensional element of line with two nodes 
(1) Solutions: Omission. 
(2) Using the equations of the shape functions, compute the following derivatives: 


Solutions: 
The shape functions: 


zt 1 — - 1+€ 
mo=a@=E 2-H moa - =! 
The derivatives: 
ON 1 ON, 1 ON | PN ud 
OE LVE L? BB” 9B 


4.6.2 Two-dimensional element of triangle with three nodes. 


(1) Solutions: Omission. 
(2) Using the equations of the shape functions, compute the following derivatives: 
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Solutions: 
The shape functions: 


1 
N(x, y) = 5A (at baz + cay), a=1,2,3 


The derivatives: 


ON, b, ON, Ca UN, QN, g Na 
Ox 2A’ Oy 24° ðry "Oy ^" 


4.6.3 'T'wo-dimensional element of rectangle with four nodes 


(1) Solutions: Omission. 
(2) Using the equations of the shape functions, compute the following derivatives: 


Solutions: 
The shape functions: 


1 1 
N,-40- 890-7, Ny =7(1+6)(1 —n) 
1 1 
N =30 480+ n), N =30 - 6)(1 +) 
The derivatives: 
OM, 1 ON, — ON, 1 ON 1 l 
UN do a e pa SONG ea a 


ON, 1 EM 1 ON, 1 N, 1 
dén 4’ O&m 4’ O&n 4’ 0Oéq 4 


ON QN QN QN 


1 , ki 0 
ae ae ae ae 


SN QN QN QN 


0 
On? "Og? ” an? "Og 


4.6.4 Three-dimensional element of tetrahedron with four nodes. 


(1) Solutions: Omission. 
(2) Using the equations of the shape functions, compute the following derivatives: 


Solutions: 
'The shape function: 
1 


gy (ae + bat + cayt daz), a=1, 2, 3,4 


N,(2, Y, z) = 
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The derivatives: 


ON, b ON, Ca ON, d, ON, , ON, , ON, 


Ox  6V' Oy 6V’ Oz 06V ary — 07 Oyz pis Ozu ni 
ON, Q9N, 9 EM ON o 
Oxyz | Ott "Oy "o2 


4.6.5 Three-dimensional element of hexahedron with eight nodes. 
(1) Solutions: Omission. 
(2) Using the equations of the shape functions, compute the following derivatives: 


Solutions: 
The shape function: 


1 
NS, n, 0 = z (1+ aS) (1+ num (1+ Caf) 
The derivatives: 

ON, 1 ON, 

dé uri (1 E 5,1) (12 CC); ôn =" Nall + 6, ENL A455, 
ON, a! PN, i 

at NM ball + ése) NaN), sn (1 1 CaS) 
Pl Ede. SG + nat) 

ant 8" rae i 
ON, >` ON, = ON, A Na 1 

ae = 0, On? a 0, ad — 0, dén C CE Caleed 


4.6.6 Summarise the functions of Lagrange interpolation for constructing one-, two-, 
and three-dimensional elements. 
Solutions: Omission. 


Chapter 5 


5.3.1 Based on the parametric mapping defined by the relation between the coor- 
dinate interpolation and the dependent-variable interpolation for a one-dimensional 
Lagrange element, describe sub-parametric, isoparametric, and super-parametric 
interpolation. 

Solutions: 

The types of parametric mapping that can be considered in an analysis are defined 
by the relation between the coordinate interpolation and the dependent-variable 
interpolation. We denote the mapping by: 


aë = X NAE) z = Ni(€) ap +.No(€) ay+---+Nn(€) ef, -1SES1 — (52a) 
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a= X N,(D à; = N (E) f+ M(E) O+---+Np(6) H%,-1<6<1 — (5.2b) 


a=1 


Once the shape functions are available in terms of the parent coordinates, we may 
immediately use the concept of parametric mapping. We have the following three 
cases: 


(1) Sub-parametric interpolation: The order n of the interpolation for z is lower 
than that order m for u. 

(2) Isoparametric interpolation: The order n of the interpolation for x is the same 
as that order m for u. 

(3) Super-parametric interpolation: The order n of the interpolation for z is higher 
than that order m for u. 


9.3.2 Consider the two-dimensional rectangle isoparametric element as shown in 
figure 5.6. 


(a) Write the expression for an isoparametric mapping of coordinates in this 
element. 

(b) Determine the location of the local coordinates č and s, which define the cen- 
troid of this element. 

(c) Compute the expression for the Jacobian matrix J of this element, and evaluate 
the Jacobian at the centroid. 

(d) Compute the derivatives of the shape function N; at the centroid. 


Solutions: 


(a) Write the expression for an isoparametric mapping of coordinates in this 
element. 
'The shape functions are 


Me eden, Avenged 


Ns = 0+8 +n), M=- En) 


aim ALE 
mle Ale 


The mapping between the global and the local coordinate system is expressed 
as: 


a(é,n) = >> Neam. (1+ ze us 
ü—90-—1. 54 (0 -90-m.,, (14- 614-1) 
4 ] 4 T 4 
P m 


x 3 


0 


3(3 — n) +é) 
4 
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En = Y Na) = YO eS te, 


(Lo) (ba) x04 tO) xta (1+ &)(1 1) " 
4 4 4 
,U-90*»,, 
3(1 +n) 
2 


3 


(b) Determine the location of the local coordinates € and y, which define the cen- 
troid of this element. 
The centre of the irregular quadrilateral is located at: 


_ 0+6+3+0 9 
ME 4 ~4 
— 0+0+3+3 3 
y= = 


(+7) 


(c) Compute the expression for the Jacobian matrix J of this element and evaluate 
the Jacobian matrix at the centroid. 
For this two-dimensional problem, the Jacobian is. 


Oxdy 9 1 
ORE | |a (: -3 n) : 
Oxdy 3 
anon cu 9 


The Jacobian at the centroid is 


NI% =) 
w 
NIW © 
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(d) Compute the derivatives of the shape function N3 at the centroid. 


The derivatives of the shape function N3 at the centroid are 


ON. ON; cr Dig 4 qi 1 
ðs | _ i) 96. [9 4 | {9 AL. ]9 
aN[ ƏN | |2 -pree[ qe c3ppae[f 2 
Oy an 9 3 4 9 3144 9 


using the inverse of the Jacobian matrix 


3 4 

=. d ME V 
ji: 8]? | |9 

2713 9 2 2 

4 4 9 3 


5.3.3 Compute stiffness X43 for one-dimensional element with three nodes using the 
following manners 


(a) Numerical integration by one, two, and three integration points, respectively. 
(b) Direct integration. 


Solutions: 
(a) Numerical integration 
'The shape functions and Jacobian of one-dimensional element with three nodes 


are 


1 1 
M=Z6E-1), Ma D Male 


1 1 1 
TE (e-z) (e+ 5) af tin = jet Sal o) 


For the case where the coordinate for node 3 is centered between nodes 1 and 2, 
the derivatives are linear in £, and the Jacobian is constant Ih. 
The stiffness for the three-node element involves the integral of 


pe A aN, s, [5 aN, a i 
0 Or Ox Ox 


[63 C363 as 
eae (C3 (yeu: 


(e-i)  (ei)e:9 — ew 
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Therefore, the stiffness parameter to be obtained is 


T 
a= f (e- 5) (2948 


Using a one-point quadrature formula, we can obtain that 


1 
Kaa f (s aK 254t =F (0-5)(-2x0) x2=0 


Using a two-point quadrature formula, we can obtain that 


ka= f (: 2 ;) (—28)ae 
A IB)9Q)* 


8E, 
3h, 


Using a three-point quadrature formula, we can obtain that 


1 
Ka =f (c7 5) ooa 


-E (s) Co (ni e oao 


(b) Direct integration 


2E, f! 1 8E, 
Kaa f. (t- 5) cont = - 5; 


9.3.4 For two-dimensional rectangle element with four nodes by numerical integra- 
tion and direct integration, respectively, compute the following integration 


i ve (+ n)dčdn 
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Solutions: 
(1) Numerical integration 


1 1 
I, i (C+ n)dčdn = = ` i (Em + Nn) Wm Wn 


m=1 n=1 


= (& m )wy wy + (6i +g) wi we + (So + m; )ws wi + (E2 + N2) wo we 
epee gegen) 
xixi+ (G+ Aun 


=0 


(2) Direct integration 


1 


f. [imm f Een] a f 2nd = 0 


Chapter 6 


6.5.1 Derive the weak form for general elasticity problems by following the steps 
below: 


(a) Multiply each equation by an appropriate arbitrary function. 

(b) Integrate this product over the space domain of the problem. 

(c) Use integration by parts to reduce the order of the derivatives to a minimum. 
(d) Introduce boundary conditions if possible. 


Solutions: Omission. 


6.5.2 In example 6.1, compute the value of the last coefficient in the stiffness matrix 
K}, of element 1: 


E [ON3 ƏN; ON3 ON3 
1 1 —2v)/2 Q 
HP ay | ay Ox s | ofa ba 
Solutions: 


The derivatives of shape functions (as shown in example 6.1) are 


ON; 1 

Or | J 30 Urn) 
ON; { -— . 
Oy 250 * 9 


75 
and the determinant of the Jacobian matrix is j BE By substituting these 


formulas into the integral equation, we can get 
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ieee ud p i 2| fan 


-f asas (35) care -9+ (4) a«wa- ae jačan 


ü L a 2 x 0.3) (35) eo x0.7 + (=) [IL P^ atas 


= 1923076.92 i 


Qe 


[0.065625 - (1 + é)? + 0.008333 - (1 + n] dédn 


O 1\? 
0.065625 «| 1 — = | +0.008333- |1- — ] |-1-1+ 


1\2 
0.065625 - | 1— =] + 0.008333 - (1+ 5) |: SEE 


E 
= 1923076.92 - 2 
1 1 
0.065625 - (1 ) 0.008333 - (1- a .-1-1+ 
v3 V3 
i? iM 
0.065625 - (: + 3 + 0.008333 - (1+ 45) | -1 
= 758543.6 
Chapter 7 
7.2.1 Using the LU decomposition of the stiffness matrix K, compute the 
displacements: 
6 2 1 8 
ku =f with K=}2 6 2 f^44 
126 1 
Solutions: 


Compute the coefficients in the matrix in the forward elimination process: 


Uj—-6 Ug=2 Ug-l 


1 16 5 1 5 
L = — — É— L n L = — 9 = — 
a= 3 Us» 3 U23 p hisp MSp U33 16 


The lower triangular matrix L and upper triangular matrix U are 


1 0 0 6 2 1 
— 1/3 1 o|],U=|0 16/3 5/3 
1/6 5/16 1 0 0 85/16 
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Compute the displacements solutions in the back substitution process: 
Ly=f and Uu=y 
The components of solutions y are 
u-h-8 
2-1 
1 4 
B-h- la h- Ian =4-73X8=7 


ca 5 4 3 


V 1 
V = f 3 3j Yj f 3101 32 12 6 1653 1 


So the following equation holds true: 


vù =y= [s 


The components of solutions u are 


Y3 3 T 16 12 


Uu E 
bf. 2*8. 1 wD) = 107 
6 17 85 ~ 85 


The solutions u are 


a= [E 5 - 
85 17 85 


The above solutions are correct by checking 


107 
85 
pem ET 8 
Ká—-|2 6 2 EN 
i Deel} E 1 
12 
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Chapter 8 


8.3.1 Summarise the basic methods and procedures in adaptive analysis. 
Solutions: 


(1) Aversion refinement: The same order of element polynomial is used, but the 
elements are changed in size (in some locations, they are made larger, and in 
others, smaller) to achieve maximum efficiency in reaching the desired solution. 

(2) p-version refinement: The procedure uses the same element size and simply 
increases the order of the element polynomials. 

(3) hp-version refinement: The procedure combines the h- and p-version refine- 
ments. In this procedure, both the element size and the degree of the element 
polynomial are altered. 


Chapter 9 


9.4.1 Using the one-, two-, and three-dimensional programs provided above, com- 
pute the displacements and stresses for the following one-, two-, and 
three-dimensional models with various values of the parameters a, b, and c, as shown 
in figure 9.51. 

Solutions: Omission. 


9.4.2 Using the analytical results of the theoretical models, verify the correctness of 
the programs in this chapter, and discuss the influence of different mesh densities on 
the accuracy of the solutions. 

Solutions: Omission. 


9.4.3 Based on the one-dimensional program of beam deformation, use shape 
functions, and numerical integration to compute the results. 
Solutions: Omission. 


9.4.4 Based on the two-dimensional program of plane strain problem, compute the 
solutions for the plane stress model. 
Solutions: 

The difference between the governing equations of plane stress and plane strain 
problems lies in the constitutive equation. If the physical parameters in the con- 
stitutive equation are transformed as follows, the two problems are equivalent: 


E— E(1+2v)/(1+v) 


From plane strain to plan stress: 
y — v/(14- v) 


E > E/üG-r 


From plane stress to plan strain: E 
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Therefore, by changing the above physical parameters in the program, the 
computation and analysis of the plane stress model can be implemented according to 
the plane strain program. 
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